Modify (d/z)rand(s/a) and the tests

5 files changed, 129 insertions, 44 deletions

diff --git a/src/auxiliaryFunctions/rand/drands.c b/src/auxiliaryFunctions/rand/drands.c
index 5ce70945..9b7c4cf4 100644
--- a/src/auxiliaryFunctions/rand/drands.c
+++ b/src/auxiliaryFunctions/rand/drands.c
@@ -1,6 +1,6 @@
 /*
  *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- *  Copyright (C) 2007-2008 - INRIA - Bruno JOFRET
+ *  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
@@ -9,11 +9,63 @@
  *  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) {
-  /* FIXME : Dummy ... */
-  return 0;
+	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;
 }
 
+
diff --git a/src/auxiliaryFunctions/rand/testRand.c b/src/auxiliaryFunctions/rand/testRand.c
index f231c95a..c700df8e 100644
--- a/src/auxiliaryFunctions/rand/testRand.c
+++ b/src/auxiliaryFunctions/rand/testRand.c
@@ -12,68 +12,84 @@
 
 #include "testRand.h"
 
-int srandsTest() {
-
-  float nan = 0.0f / 0.0f;
-  float result = nan;
-
-  printf(">> Float \n");
-  result = srands();
-
-  assert(&result != &nan);
-  return 0;
-}
+/*
 
 int drandsTest() {
 
   double nan = 0.0 / 0.0;
   double result = nan;
-
+  int i=0;
   printf(">> Double\n");
+  for (i=0;i<10;i++){
   result = drands();
-
+  printf("%1.20f\n",result);
+  }
+	
   assert(&result != &nan);
 
   return 0;
 }
+*/
 
-int crandsTest() {
-
-  floatComplex nan_nan = FloatComplex(0.0f/0.0f, 0.0f/0.0f);
-  floatComplex result = nan_nan;
+/*
+int zrandsTest() {
 
-  printf(">> Float Complex\n");
-  result = crands();
+  double nan = 0.0 / 0.0;
+  doubleComplex result = DoubleComplex(nan,nan);
+  int i=0;
+  printf(">> DoubleComplex\n");
+  for (i=0;i<10;i++){
+  result = zrands();
+  printf("%1.20f+%1.20f\n",zreals(result),zimags(result));
+  }
 
-  assert(&result != &nan_nan);
+  assert(zreals(result) != nan);
 
   return 0;
 }
+*/
 
-int zrandsTest() {
-  doubleComplex nan_nan = DoubleComplex(0./0., 0./0.);
-  doubleComplex result = nan_nan;
 
-  printf(">> Double Complex\n");
-  result = zrands();
+int drandaTest() {
+  int i;
+  double nan = 0.0 / 0.0;
+  double *result;
+  result=malloc((uint)12*sizeof(double));
+  printf(">> Double Array\n");
+  dranda(result,12);
+  for (i=0;i<12;i++){
+  printf("%1.20f\n",result[i]);
 
-  assert(&result != &nan_nan);
+  assert(result[i] != nan);}
 
   return 0;
 }
+  
+
+/*
+int zrandaTest() {
+  int i;
+  double nan = 0.0 / 0.0;
+  doubleComplex *result;
+  result=malloc((uint)12*sizeof(doubleComplex));
+  printf(">> DoubleComplex Array\n");
+  zranda(result,12);
+  for (i=0;i<12;i++){
+  printf("%1.20f+%1.20f\n",zreals(result[i]),zimags(result[i]));
+
+  assert(zreals(result[i]) != nan);}
 
+  return 0;
+}
+*/
 
 int testRand() {
-  int srandsTestStatus, drandsTestStatus = 0;
-  int crandsTestStatus, zrandsTestStatus = 0;
+/*  int drandsTestStatus = 0;*/
+  int zrandaTestStatus = 0;
   printf("\n>>>> Rand Tests\n");
-  srandsTestStatus = srandsTest();
-  drandsTestStatus = drandsTest();
-  crandsTestStatus = crandsTest();
-  zrandsTestStatus = zrandsTest();
-
-  return (srandsTestStatus + drandsTestStatus +
-	  crandsTestStatus + zrandsTestStatus);
+ /* drandsTestStatus = drandsTest();*/
+  zrandaTestStatus = drandaTest();
+  return (zrandaTestStatus);
 }
 
 int main(void) {
diff --git a/src/auxiliaryFunctions/rand/testRand.h b/src/auxiliaryFunctions/rand/testRand.h
index b55a0162..48f05788 100644
--- a/src/auxiliaryFunctions/rand/testRand.h
+++ b/src/auxiliaryFunctions/rand/testRand.h
@@ -22,11 +22,11 @@
 int srandsTest(void);
 
 int drandsTest(void);
-
+int drandaTest(void);
 int crandsTest(void);
 
 int zrandsTest(void);
-
+int zrandaTest(void);
 int testRand(void);
 
 #endif /* ! __TESTRAND_H__ */
diff --git a/src/auxiliaryFunctions/rand/zranda.c b/src/auxiliaryFunctions/rand/zranda.c
index 7dd04a16..98cc10cf 100644
--- a/src/auxiliaryFunctions/rand/zranda.c
+++ b/src/auxiliaryFunctions/rand/zranda.c
@@ -11,10 +11,23 @@
  */
 
 #include "rand.h"
-
+#include <stdio.h>
 void zranda(doubleComplex *out, int size) {
   int i = 0;
+  double *Real_randa, *Imag_randa;
+  Real_randa = malloc((uint)size*sizeof(double));
+  Imag_randa = malloc((uint)size*sizeof(double));
+   
   for (i = 0 ; i < size ; ++i) {
-    out[i] = zrands();
+    Real_randa[i]=drands();
+  }
+
+  for (i = 0 ; i < size ; ++i) {  
+    Imag_randa[i]=drands();
   }
+
+  for (i = 0 ; i < size ; ++i) {  
+    out[i]=DoubleComplex(Real_randa[i],Imag_randa[i]);
+  }	
+
 }
diff --git a/src/auxiliaryFunctions/rand/zrands.c b/src/auxiliaryFunctions/rand/zrands.c
index e9b2bd33..3011545f 100644
--- a/src/auxiliaryFunctions/rand/zrands.c
+++ b/src/auxiliaryFunctions/rand/zrands.c
@@ -13,6 +13,10 @@
 #include "rand.h"
 
 doubleComplex zrands(void) {
-  /* FIXME : Implementation */
-  return DoubleComplex(0, 0);
+  double Real_rand=drands();
+  double Imag_rand=drands();
+  return DoubleComplex(Real_rand,Imag_rand );
+  
+  /* We can't do return DoubleComplex(drands(),drands()) 
+  cause the result is false; it give imag+real*i instead of real+imag*i */
 }