1 files changed, 71 insertions, 0 deletions
diff --git a/2.3-1/src/c/auxiliaryFunctions/rand/drands.c b/2.3-1/src/c/auxiliaryFunctions/rand/drands.c
new file mode 100644
index 00000000..9b7c4cf4
--- /dev/null
+++ b/2.3-1/src/c/auxiliaryFunctions/rand/drands.c
@@ -0,0 +1,71 @@
+/*
+ *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ *  Copyright (C) 2008 - INRIA - Arnaud TORSET
+ *
+ *  This file must be used under the terms of the CeCILL.
+ *  This source file is licensed as described in the file COPYING, which
+ *  you should have received as part of this distribution.  The terms
+ *  are also available at
+ *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ *
+ */
+ 
+ 
+/*
+   URAND, A UNIVERSAL RANDOM NUMBER GENERATOR
+   BY, MICHAEL A. MALCOLM, CLEVE B. MOLER,
+   STAN-CS-73-334, JANUARY 1973,
+   COMPUTER SCIENCE  DEPARTMENT,
+   School of Humanities and Sciences, STANFORD UNIVERSITY,
+   ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/73/334/CS-TR-73-334.pdf
+
+*/
+
+
+
+#include "rand.h"
+#include <stdio.h>
+
+double	drands(void) {
+	int m=1;
+	const int itwo=2;
+	static int m2=0,halfm=0,ia=0, ic=0, mic=0, iy=0;
+	static double s=0.0;
+  
+      if (m2==0){
+	/*  if first entry, compute machine integer word length */
+		while (m>m2){
+		m2=m;
+		m=itwo*m2;
+		}
+		halfm = m2;
+	
+		/* compute multiplier and increment for linear congruential method */
+		ia = 8*(int)(halfm*atan(1.0)/8.0) + 5;
+		ic = 2*(int)(halfm*(0.5-sqrt(3.0)/6.0)) + 1;
+		mic = (m2 - ic) + m2;
+	
+		/* s is the scale factor for converting to floating point */
+		s = 0.5/halfm;
+	}
+	
+	/* compute next random number	*/
+	iy = iy*ia;
+
+	/* the following statement is for computers which do not allow
+	   integer overflow on addition  */
+      if (iy > mic) iy = (iy - m2) - m2;
+
+      iy = iy + ic;
+
+	/* the following statement is for computers where the
+	   word length for addition is greater than for multiplication */
+      if (iy/2 > m2) iy = (iy - m2) - m2;
+
+	/* the following statement is for computers where integer
+	   overflow affects the sign bit */
+      if (iy < 0) iy = (iy + m2) + m2;
+      return (double)iy*s;
+}
+
+