latest version

1 files changed, 153 insertions, 0 deletions

diff --git a/graph-checker/nauty/naurng.c b/graph-checker/nauty/naurng.c
new file mode 100644
index 0000000..dcd4d2c
--- /dev/null
+++ b/graph-checker/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()); 
+}