redo random number generation

wallace.c, main.c, randnumb.c, sharedspice.c:
set a seed value by getpid, may be overridden by 'set rndseed=val'
add func gauss1() returning only a single value
speed up gauss0()
add checkseed() to rgauss() and poisson() to read 'rndseed', if changed
This commit is contained in:
h_vogt 2016-08-04 19:27:28 +02:00 committed by rlar
parent 802c7c3261
commit 5aabd5e8f2
5 changed files with 56 additions and 39 deletions

View File

@ -79,10 +79,12 @@ initw(void)
double totsqr, nomsqr;
unsigned int coa;
#ifdef HasMain
/* initialize the uniform generator */
srand((unsigned int) getpid());
// srand(17);
TausSeed();
#endif
ScaleGauss = 1.;
newpools = 1;

View File

@ -7,6 +7,7 @@ extern void initw(void);
extern void checkseed(void); /* seed random or set by 'set rndseed=value'*/
extern double drand(void);
extern double gauss0(void);
extern double gauss1(void);
extern int poisson(double);
extern double exprand(double);

View File

@ -887,8 +887,9 @@ main(int argc, char **argv)
cp_program = ft_sim->simulator;
srand((unsigned int) getpid());
TausSeed();
int ii = getpid();
cp_vset("rndseed", CP_NUM, &ii);
checkseed();
/* --- Process command line options --- */
for (;;) {
@ -1188,15 +1189,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)) {
time_t acttime = time(NULL);
rseed = (unsigned int) acttime;
}
srand(rseed);
initw();
}
initw();
#endif
if (!ft_servermode) {

View File

@ -2,8 +2,11 @@
Copyright 2008 Holger Vogt
**********/
/* Care about random numbers
The seed value is set as random number in main.c, line 746.
A fixed seed value may be set by 'set rndseed=value'.
The seed value is obtained from getpid() in main.c upon start-up
and stored in a variable named 'rndseed'.
A fixed seed value may be set by 'set rndseed=value' in .spiceinit
or interactively. This value then is acknowledged by the next call
to checkseed().
*/
@ -68,9 +71,13 @@ void rgauss(double* py1, double* py2);
/* 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 getpid will be used, if no 'set rndseed=val' is given.
Checkseed should be put in front of any call to rand or CombLCGTaus
to catch any change of rndseed.
Checkseed may not be used to reset the random number generator, if
rndseed has not been changed.
*/
void checkseed(void)
{
@ -81,10 +88,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); */
}
}
@ -190,16 +197,18 @@ unsigned int CombLCGTausInt2(void)
}
/*** gauss ***/
/*** gauss ***
for speed reasons get two values per pass */
double gauss0(void)
{
static bool gliset = TRUE;
static double glgset = 0.0;
double fac,r,v1,v2;
if (gliset) {
checkseed();
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); */
@ -213,6 +222,24 @@ double gauss0(void)
}
}
/*** gauss ***
to be reproducible, we just use one value per pass */
double gauss1(void)
{
double fac, r, v1, v2;
checkseed();
do {
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); */
fac = sqrt(-2.0 * log(r) / r);
return v2 * fac;
}
/* Polar form of the Box-Muller generator for Gaussian distributed
random variates.
Generator will be fed with two uniformly distributed random variates.
@ -222,14 +249,14 @@ double gauss0(void)
void rgauss(double* py1, double* py2)
{
double x1, x2, w;
checkseed();
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;
@ -243,6 +270,7 @@ int poisson(double lambda)
{
int k=0; //Counter
const int max_k = 1000; //k upper limit
checkseed();
double p = CombLCGTaus(); //uniform random number
double P = exp(-lambda); //probability
double sum=P; //cumulant

View File

@ -582,8 +582,9 @@ ngSpice_Init(SendChar* printfcn, SendStat* statusfcn, ControlledExit* ngspiceexi
/* program name*/
cp_program = ft_sim->simulator;
srand((unsigned int) getpid());
TausSeed();
int ii = getpid();
cp_vset("rndseed", CP_NUM, &ii);
checkseed();
/*parameter fetcher, used in show, alter, altermod */
if_getparam = spif_getparam_special;
@ -648,15 +649,7 @@ bot:
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)) {
time_t acttime = time(NULL);
rseed = (unsigned int) acttime;
}
srand(rseed);
initw();
}
initw();
#endif
// com_version(NULL);