From 2786fcb86eef3280915b361c50ae5d348b950d76 Mon Sep 17 00:00:00 2001 From: Holger Vogt Date: Mon, 2 Jul 2018 21:18:28 +0200 Subject: [PATCH] New .option seed=[val|random] --- src/frontend/numparam/xpressn.c | 4 +-- src/main.c | 10 +------- src/maths/cmaths/cmath2.c | 2 +- src/maths/misc/randnumb.c | 43 +++++++++++++++------------------ 4 files changed, 24 insertions(+), 35 deletions(-) diff --git a/src/frontend/numparam/xpressn.c b/src/frontend/numparam/xpressn.c index 40c4eaad5..ed957fd7b 100644 --- a/src/frontend/numparam/xpressn.c +++ b/src/frontend/numparam/xpressn.c @@ -47,7 +47,7 @@ agauss(double nominal_val, double abs_variation, double sigma) { double stdvar; stdvar = abs_variation / sigma; - return (nominal_val + stdvar * gauss0()); + return (nominal_val + stdvar * gauss1()); } @@ -56,7 +56,7 @@ gauss(double nominal_val, double rel_variation, double sigma) { double stdvar; stdvar = nominal_val * rel_variation / sigma; - return (nominal_val + stdvar * gauss0()); + return (nominal_val + stdvar * gauss1()); } diff --git a/src/main.c b/src/main.c index 55e1c931b..2b213eb3b 100644 --- a/src/main.c +++ b/src/main.c @@ -1177,15 +1177,7 @@ main(int argc, char **argv) fprintf(cp_out, "SoS %f, seed value: %ld\n", renormalize(), rseed); } #elif defined(WaGauss) - { - unsigned int rseed = 66; - if (!cp_getvar("rndseed", CP_NUM, &rseed, 0)) { - time_t acttime = time(NULL); - rseed = (unsigned int) acttime; - } - srand(rseed); - initw(); - } + initw(); #endif if (!ft_servermode) { diff --git a/src/maths/cmaths/cmath2.c b/src/maths/cmaths/cmath2.c index 85a1ea8f0..4b35335df 100644 --- a/src/maths/cmaths/cmath2.c +++ b/src/maths/cmaths/cmath2.c @@ -305,7 +305,7 @@ cx_sgauss(void *data, short int type, int length, int *newlength, short int *new d = alloc_d(length); *newtype = VF_REAL; for (i = 0; i < length; i++) { - d[i] = gauss0(); + d[i] = gauss1(); } return ((void *) d); } diff --git a/src/maths/misc/randnumb.c b/src/maths/misc/randnumb.c index b033b92d2..008598273 100644 --- a/src/maths/misc/randnumb.c +++ b/src/maths/misc/randnumb.c @@ -69,9 +69,10 @@ static bool seedinfo = FALSE; /* Check if a seed has been set by the command 'set rndseed=value' - in spinit with integer value > 0. If available, call srand(value). - This will override the call to srand in main.c. - Checkseed should be put in front of any call to rand or CombLCGTaus.. . + in spinit, .spiceinit or in a .control section + with integer value > 0. If available, call srand(value). + rndseed set in main.c to 1, if no 'set rndseed=val' is given. + Called from functions in cmath2.c. */ void checkseed(void) { @@ -82,10 +83,10 @@ void checkseed(void) if ((newseed > 0) && (oldseed != newseed)) { srand((unsigned int)newseed); TausSeed(); + if (oldseed > 0) /* no printout upon start-up */ + printf("Seed value for random number generator is set to %d\n", newseed); oldseed = newseed; - printf("Seed value for random number generator is set to %d\n", newseed); } -/* else printf("Oldseed %d, newseed %d\n", oldseed, newseed); */ } } @@ -93,8 +94,6 @@ void checkseed(void) /* uniform random number generator, interval [-1 .. +1[ */ double drand(void) { - checkseed(); -// return ( 2.0*((double) (RR_MAX-abs(rand())) / (double)RR_MAX-0.5)); return 2.0 * CombLCGTaus() - 1.0; } @@ -191,8 +190,8 @@ unsigned int CombLCGTausInt2(void) } -/*** gauss ***/ - +/*** gauss *** + for speed reasons get two values per pass */ double gauss0(void) { static bool gliset = TRUE; @@ -200,7 +199,8 @@ double gauss0(void) double fac,r,v1,v2; if (gliset) { do { - v1 = drand(); v2 = drand(); + v1 = 2.0 * CombLCGTaus() - 1.0; + v2 = 2.0 * CombLCGTaus() - 1.0; r = v1*v1 + v2*v2; } while (r >= 1.0); /* printf("v1 %f, v2 %f\n", v1, v2); */ @@ -239,19 +239,19 @@ double gauss1(void) void rgauss(double* py1, double* py2) { - double x1, x2, w; + double x1, x2, w; - do { - x1 = 2.0 * CombLCGTaus() - 1.0; - x2 = 2.0 * CombLCGTaus() - 1.0; - w = x1 * x1 + x2 * x2; - } while ( w >= 1.0 ); + do { + x1 = 2.0 * CombLCGTaus() - 1.0; + x2 = 2.0 * CombLCGTaus() - 1.0; + w = x1 * x1 + x2 * x2; + } while ( w >= 1.0 ); - w = sqrt( (-2.0 * log( w ) ) / w ); + w = sqrt( (-2.0 * log( w ) ) / w ); - *py1 = x1 * w; - *py2 = x2 * w; -} + *py1 = x1 * w; + *py2 = x2 * w; +} @@ -278,13 +278,11 @@ int poisson(double lambda) double exprand(double mean) { double expval; - checkseed(); expval = -log(CombLCGTaus()) * mean; return expval; } - /* seed random number generators immediately * command "setseed" * take value of variable rndseed as seed @@ -329,4 +327,3 @@ setseedinfo(void) { seedinfo = TRUE; } -