summaryrefslogtreecommitdiff
path: root/src/c/auxiliaryFunctions/rand
diff options
context:
space:
mode:
Diffstat (limited to 'src/c/auxiliaryFunctions/rand')
-rw-r--r--src/c/auxiliaryFunctions/rand/dranda.c2
-rw-r--r--src/c/auxiliaryFunctions/rand/drands.c43
2 files changed, 34 insertions, 11 deletions
diff --git a/src/c/auxiliaryFunctions/rand/dranda.c b/src/c/auxiliaryFunctions/rand/dranda.c
index f467827..3defc26 100644
--- a/src/c/auxiliaryFunctions/rand/dranda.c
+++ b/src/c/auxiliaryFunctions/rand/dranda.c
@@ -16,5 +16,5 @@ void dranda(double *out, int size) {
int i = 0;
for (i = 0 ; i < size ; ++i) {
out[i] = drands();
- }
+ }
}
diff --git a/src/c/auxiliaryFunctions/rand/drands.c b/src/c/auxiliaryFunctions/rand/drands.c
index 9b7c4cf..658cd0b 100644
--- a/src/c/auxiliaryFunctions/rand/drands.c
+++ b/src/c/auxiliaryFunctions/rand/drands.c
@@ -34,37 +34,60 @@ double drands(void) {
if (m2==0){
/* if first entry, compute machine integer word length */
- while (m>m2){
+ while (m>m2){
m2=m;
m=itwo*m2;
+
}
+
halfm = m2;
-
- /* compute multiplier and increment for linear congruential method */
+
+ /* 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 */
+
+ 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;
+ 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;
+ 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;
+ if (iy < 0)
+ {
+ iy = (iy + m2) + m2;
+
+
+ }
+
+
return (double)iy*s;
}