diff options
Diffstat (limited to 'nauty/naurng.c')
| -rw-r--r-- | nauty/naurng.c | 153 |
1 files changed, 153 insertions, 0 deletions
diff --git a/nauty/naurng.c b/nauty/naurng.c new file mode 100644 index 0000000..dcd4d2c --- /dev/null +++ b/nauty/naurng.c @@ -0,0 +1,153 @@ +/* naurng.c + * + This file contains the code for a high-quality random number + generator written by Don Knuth. The auxilliary routine + ran_arr_cycle() has been modified slightly, and ran_init() is new. + + To use it: + + 0. #include "naurng.h" (or "naututil.h" if you are using nauty) + + 1. Call ran_init(seed), where seed is any long integer. + This step is optional, but if you don't use it you + will always get the same sequence of random numbers. + + 2. For each random number, use the NEXTRAN macro. It will + give a random value in the range 0..2^30-1. Alternatively, + KRAN(k) will have a random value in the range 0..k-1. + KRAN(k) actually gives you NEXTRAN mod k, so it is not + totally uniform if k is very large. In that case, you + can use the slightly slower GETKRAN(k,var) to set the + variable var to a better random number from 0..k-1. + + Brendan McKay, July 2002. Fixed Nov 2002 on advice of DEK. + Nov 2022. Added ran_init_time(). +*/ + +/* This program by D E Knuth is in the public domain and freely copyable + * AS LONG AS YOU MAKE ABSOLUTELY NO CHANGES! + * It is explained in Seminumerical Algorithms, 3rd edition, Section 3.6 + * (or in the errata to the 2nd edition --- see + * http://www-cs-faculty.stanford.edu/~knuth/taocp.html + * in the changes to Volume 2 on pages 171 and following). */ + +/* N.B. The MODIFICATIONS introduced in the 9th printing (2002) are + included here; there's no backwards compatibility with the original. */ + +/* If you find any bugs, please report them immediately to + * taocp@cs.stanford.edu + * (and you will be rewarded if the bug is genuine). Thanks! */ + +/************ see the book for explanations and caveats! *******************/ +/************ in particular, you need two's complement arithmetic **********/ + +#include "naurng.h" + +#define KK 100 /* the long lag */ +#define LL 37 /* the short lag */ +#define MM (1L<<30) /* the modulus */ +#define mod_diff(x,y) (((x)-(y))&(MM-1)) /* subtraction mod MM */ + +static TLS_ATTR long ran_x[KK]; /* the generator state */ + +static void +ran_array(long aa[],int n) +{ + int i,j; + for (j=0;j<KK;j++) aa[j]=ran_x[j]; + for (;j<n;j++) aa[j]=mod_diff(aa[j-KK],aa[j-LL]); + for (i=0;i<LL;i++,j++) ran_x[i]=mod_diff(aa[j-KK],aa[j-LL]); + for (;i<KK;i++,j++) ran_x[i]=mod_diff(aa[j-KK],ran_x[i-LL]); +} + +/* the following routines are from exercise 3.6--15 */ +/* after calling ran_start, get new randoms by, e.g., "x=ran_arr_next()" */ + +#define RNG_QUALITY 1009 /* recommended quality level for high-res use */ +static TLS_ATTR long ran_arr_buf[RNG_QUALITY]; +static long ran_arr_dummy=-1; +static long ran_arr_started=-1; +static TLS_ATTR long *ran_arr_ptr = &ran_arr_dummy; + +#define TT 70 /* guaranteed separation between streams */ +#define is_odd(x) ((x)&1) /* units bit of x */ + +static void +ran_start(long seed) +{ + int t,j; + long x[KK+KK-1]; /* the preparation buffer */ + long ss=(seed+2)&(MM-2); + + for (j=0;j<KK;j++) { + x[j]=ss; /* bootstrap the buffer */ + ss<<=1; if (ss>=MM) ss-=MM-2; /* cyclic shift 29 bits */ + } + x[1]++; /* make x[1] (and only x[1]) odd */ + for (ss=seed&(MM-1),t=TT-1; t; ) { + for (j=KK-1;j>0;j--) x[j+j]=x[j], x[j+j-1]=0; /* "square" */ + for (j=KK+KK-2;j>=KK;j--) + x[j-(KK-LL)]=mod_diff(x[j-(KK-LL)],x[j]), + x[j-KK]=mod_diff(x[j-KK],x[j]); + if (is_odd(ss)) { /* "multiply by z" */ + for (j=KK;j>0;j--) x[j]=x[j-1]; + x[0]=x[KK]; /* shift the buffer cyclically */ + x[LL]=mod_diff(x[LL],x[KK]); + } + if (ss) ss>>=1; else t--; + } + for (j=0;j<LL;j++) ran_x[j+KK-LL]=x[j]; + for (;j<KK;j++) ran_x[j-LL]=x[j]; + for (j=0;j<10;j++) ran_array(x,KK+KK-1); /* warm things up */ + ran_arr_ptr=&ran_arr_started; +} + +void +ran_init(long seed) +/* Added by BDM: use instead of ran_start. */ +{ + ran_start((unsigned long)seed % (MM-2)); +} + +long +ran_init_time(long extra) +/* Added by BDM: use the real time and the argument to initialise. + Returns the value of the seed, so the same sequence can be + obtained again by calling ran_init(seed). +*/ +{ + double t; + nauty_counter ul; /* 64-bit unsigned */ + long seed; + REALTIMEDEFS + + t = NAUTYREALTIME; + if (t > 1660000000.0) ul = (nauty_counter)(t*2100001.0); + else ul = (nauty_counter)(t+212300021.0); + ul ^= (nauty_counter)(extra) * 997; + seed = (long)ul; + ran_init(seed); + + return seed; +} + +static long +ran_arr_cycle(void) +/* Modified by BDM to automatically initialise + if no explicit initialisation has been done */ +{ + if (ran_arr_ptr==&ran_arr_dummy) + ran_start(314159L); /* the user forgot to initialize */ + + ran_array(ran_arr_buf,RNG_QUALITY); + + ran_arr_buf[KK]=-1; + ran_arr_ptr=ran_arr_buf+1; + return ran_arr_buf[0]; +} + +long +ran_nextran(void) +{ + return (*ran_arr_ptr>=0 ? *ran_arr_ptr++ : ran_arr_cycle()); +} |