From ba539d9ef1dd6954ba2037882d8bd1cb7586b248 Mon Sep 17 00:00:00 2001 From: h_vogt Date: Sat, 28 Aug 2010 18:13:08 +0000 Subject: [PATCH] new fcn sgauss(), new rnd-no generator --- ChangeLog | 5 ++ src/frontend/parse.c | 1 + src/main.c | 1 + src/maths/cmaths/cmath2.c | 33 ++++++++ src/maths/cmaths/cmath2.h | 1 + src/maths/misc/randnumb.c | 157 +++++++++++++++++++++++++++++++++++++- 6 files changed, 194 insertions(+), 4 deletions(-) diff --git a/ChangeLog b/ChangeLog index bd3b8f298..eba1452c5 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,8 @@ +2010-08-29 Holger Vogt + * fteext.h, main.c, randnumb.c, parse.c, cmath2.c, cmath2.h: + new function sgauss(), new random number generator with very + high run length + 2010-08-19 Holger Vogt * xpressn.c: bug no. 3047884 fixed inp.c: prevent crash if .param is last line in input deck diff --git a/src/frontend/parse.c b/src/frontend/parse.c index 9a223a6b1..a8dd2e113 100644 --- a/src/frontend/parse.c +++ b/src/frontend/parse.c @@ -161,6 +161,7 @@ struct func ft_funcs[] = { { "atan", cx_atan } , { "norm", cx_norm } , { "rnd", cx_rnd } , + { "sgauss", cx_sgauss } , { "pos", cx_pos } , { "mean", cx_mean } , { "avg", cx_avg } , /* A.Roldan 03/06/05 incremental average new function */ diff --git a/src/main.c b/src/main.c index 38a0bcfb6..1d696bc2d 100644 --- a/src/main.c +++ b/src/main.c @@ -788,6 +788,7 @@ main(int argc, char **argv) cp_program = ft_sim->simulator; srandom(getpid()); + TausSeed(); /* --- Process command line options --- */ while (1) { diff --git a/src/maths/cmaths/cmath2.c b/src/maths/cmaths/cmath2.c index dff5af98d..5102b4a90 100644 --- a/src/maths/cmaths/cmath2.c +++ b/src/maths/cmaths/cmath2.c @@ -31,6 +31,7 @@ extern void srandom (unsigned int seed); #endif extern void checkseed(void); /* seed random or set by 'set rndseed=value'*/ +extern double gauss(void); /* from randnumb.c */ static double * d_tan(double *dd, int length) @@ -244,6 +245,38 @@ cx_rnd(void *data, short int type, int length, int *newlength, short int *newtyp } } +void * +cx_sgauss(void *data, short int type, int length, int *newlength, short int *newtype) +{ + *newlength = length; + checkseed(); + if (type == VF_COMPLEX) { + complex *c; + complex *cc = (complex *) data; + int i; + + c = alloc_c(length); + *newtype = VF_COMPLEX; + for (i = 0; i < length; i++) { + realpart(&c[i]) = gauss(); + imagpart(&c[i]) = gauss(); + } + return ((void *) c); + } else { + double *d; + double *dd = (double *) data; + int i; + + d = alloc_d(length); + *newtype = VF_REAL; + for (i = 0; i < length; i++) { + d[i] = gauss(); + } + return ((void *) d); + } +} + + /* Compute the avg of a vector. Created by A.M.Roldan 2005-05-21 */ diff --git a/src/maths/cmaths/cmath2.h b/src/maths/cmaths/cmath2.h index 16cdf668f..f2462b63f 100644 --- a/src/maths/cmaths/cmath2.h +++ b/src/maths/cmaths/cmath2.h @@ -12,6 +12,7 @@ void * cx_atan(void *data, short int type, int length, int *newlength, short int void * cx_norm(void *data, short int type, int length, int *newlength, short int *newtype); void * cx_uminus(void *data, short int type, int length, int *newlength, short int *newtype); void * cx_rnd(void *data, short int type, int length, int *newlength, short int *newtype); +void * cx_sgauss(void *data, short int type, int length, int *newlength, short int *newtype); void * cx_mean(void *data, short int type, int length, int *newlength, short int *newtype); void * cx_length(void *data, short int type, int length, int *newlength, short int *newtype); void * cx_vector(void *data, short int type, int length, int *newlength, short int *newtype); diff --git a/src/maths/misc/randnumb.c b/src/maths/misc/randnumb.c index e3a484f2f..8343df53b 100644 --- a/src/maths/misc/randnumb.c +++ b/src/maths/misc/randnumb.c @@ -6,9 +6,64 @@ Copyright 2008 Holger Vogt A fixed seed value may be set by 'set rndseed=value'. */ + +/* + CombLCGTaus() + Combined Tausworthe & LCG random number generator + Algorithm has been suggested in: + GPUGems 3, Addison Wesley, 2008, Chapter 37. + It combines a three component Tausworthe generator taus88 + (see P. L'Ecuyer: "Maximally equidistributed combined Tausworthe + generators", Mathematics of Computation, 1996, + www.iro.umontreal.ca/~lecuyer/myftp/papers/tausme.ps) + and a quick linear congruent generator (LCG), decribed in: + Press: "Numerical recipes in C", Cambridge, 1992, p. 284. + Generator has passed the bbattery_SmallCrush(gen) test of the + TestU01 library from Pierre L’Ecuyer and Richard Simard, + http://www.iro.umontreal.ca/~simardr/testu01/tu01.html +*/ + + +/* TausSeed creates three start values for Tausworthe state variables. + Uses rand() from , therefore values depend on the value of + seed in srand(seed). A constant seed will result in a reproducible + series of random variates. + + Calling sequence: + srand(seed); + TausSeed(); + // generate random variates randvar uniformly distributed in + // [0.0 .. 1.0[ by calls to CombLCGTaus(). + double randvar = CombLCGTaus(void); +*/ + + #include "ngspice.h" #include "cpdefs.h" #include "ftedefs.h" +#include +#include +#include +#ifdef _MSC_VER +#include +#else +#include +#endif +#include +#include // var. argumente + +/* Tausworthe state variables for double variates*/ +static unsigned CombState1 = 129, CombState2 = 130, CombState3 = 131; +static unsigned CombState4; /* LCG state variable */ + +/* Tausworthe state variables for double variates*/ +static unsigned CombState5 = 133, CombState6 = 135, CombState7 = 137; +static unsigned CombState8; /* LCG state variable */ + +unsigned TauS(unsigned *state, int C1, int C2, int C3, unsigned m); +unsigned LGCS(unsigned *state, unsigned A1, unsigned A2); +void TausSeed(void); + /* MINGW: random, srandom in libiberty.a, but not in libiberty.h */ @@ -35,6 +90,7 @@ void checkseed(void) if (cp_getvar("rndseed", CP_NUM, &newseed)) { if ((newseed > 0) && (oldseed != newseed)) { srandom(newseed); + TausSeed(); oldseed = newseed; printf("Seed value for random number generator is set to %d\n", newseed); } @@ -51,6 +107,100 @@ double drand(void) } + +void TausSeed(void) +{ + CombState1 = CombState2 = CombState3 = 0; + /* The Tausworthe initial states should be greater than 128 */ + while (CombState1 < 129) + CombState1 = rand(); + while (CombState2 < 129) + CombState2 = rand(); + while (CombState3 < 129) + CombState3 = rand(); + while (CombState4 < 129) + CombState4 = rand(); +#ifdef HVDEBUG + printf("\nTausworthe Double generator init states: %d, %d, %d, %d\n", + CombState1, CombState2, CombState3, CombState4); +#endif + CombState5 = CombState6 = CombState7 = 0; + /* The Tausworthe initial states should be greater than 128 */ + while (CombState5 < 129) + CombState5 = rand(); + while (CombState6 < 129) + CombState6 = rand(); + while (CombState7 < 129) + CombState7 = rand(); + while (CombState8 < 129) + CombState8 = rand(); +#ifdef HVDEBUG + printf("Tausworthe Integer generator init states: %d, %d, %d, %d\n", + CombState5, CombState6, CombState7, CombState8); +#endif +} + +unsigned TauS(unsigned *state, int C1, int C2, int C3, unsigned m) +{ + unsigned b = (((*state << C1) ^ *state) >> C2); + return *state = (((*state & m) << C3) ^ b); +} + +unsigned LGCS(unsigned *state, unsigned A1, unsigned A2) +{ + return *state = (A1 * *state + A2); +} + + +double CombLCGTaus() +{ + return 2.3283064365387e-10 * ( + TauS(&CombState1, 13, 19, 12, 4294967294UL) ^ + TauS(&CombState2, 2, 25, 4, 4294967288UL) ^ + TauS(&CombState3, 3, 11, 17, 4294967280UL) ^ + LGCS(&CombState4, 1664525, 1013904223UL) + ); +} + +unsigned int CombLCGTausInt() +{ + return ( + TauS(&CombState5, 13, 19, 12, 4294967294UL) ^ + TauS(&CombState6, 2, 25, 4, 4294967288UL) ^ + TauS(&CombState7, 3, 11, 17, 4294967280UL) ^ + LGCS(&CombState8, 1664525, 1013904223UL) + ); +} + + +float CombLCGTaus2() +{ + unsigned long b; + b = (((CombState1 << 13) ^ CombState1) >> 19); + CombState1 = (((CombState1 & 4294967294UL) << 12) ^ b); + b = (((CombState2 << 2) ^ CombState2) >> 25); + CombState2 = (((CombState2 & 4294967288UL) << 4) ^ b); + b = (((CombState3 << 3) ^ CombState3) >> 11); + CombState3 = (((CombState3 & 4294967280UL) << 17) ^ b); + CombState4 = (1664525 * CombState4 + 1013904223UL); + return ((CombState1 ^ CombState2 ^ CombState3 ^ CombState4) * 2.3283064365387e-10f); +} + + +unsigned int CombLCGTausInt2() +{ + unsigned long b; + b = (((CombState5 << 13) ^ CombState5) >> 19); + CombState5 = (((CombState5 & 4294967294UL) << 12) ^ b); + b = (((CombState6 << 2) ^ CombState6) >> 25); + CombState6 = (((CombState6 & 4294967288UL) << 4) ^ b); + b = (((CombState7 << 3) ^ CombState7) >> 11); + CombState7 = (((CombState7 & 4294967280UL) << 17) ^ b); + CombState8 = (1664525 * CombState8 + 1013904223UL); + return (CombState5 ^ CombState6 ^ CombState7 ^ CombState8); +} + + /*** gauss ***/ double gauss(void) @@ -60,7 +210,9 @@ double gauss(void) double fac,r,v1,v2; if (gliset) { do { - v1 = drand(); v2 = drand(); +// 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); */ @@ -73,6 +225,3 @@ double gauss(void) return glgset; } } - - -