summaryrefslogtreecommitdiff
path: root/support/graphbase/gb_flip.w
diff options
context:
space:
mode:
authorNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
committerNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
commite0c6872cf40896c7be36b11dcc744620f10adf1d (patch)
tree60335e10d2f4354b0674ec22d7b53f0f8abee672 /support/graphbase/gb_flip.w
Initial commit
Diffstat (limited to 'support/graphbase/gb_flip.w')
-rw-r--r--support/graphbase/gb_flip.w254
1 files changed, 254 insertions, 0 deletions
diff --git a/support/graphbase/gb_flip.w b/support/graphbase/gb_flip.w
new file mode 100644
index 0000000000..c24528f189
--- /dev/null
+++ b/support/graphbase/gb_flip.w
@@ -0,0 +1,254 @@
+% This file is part of the Stanford GraphBase (c) Stanford University 1992
+\def\title{GB\_\thinspace FLIP}
+@i boilerplate.w %<< legal stuff: PLEASE READ IT BEFORE MAKING ANY CHANGES!
+
+@* Introduction. This is |gb_flip|, the module used by GraphBase
+programs to generate random numbers.
+
+To use the routines in this file, first call the function |gb_init_rand(seed)|.
+Subsequent uses of the macro |gb_next_rand()| will then return pseudo-random
+integers between 0 and $2^{31}-1$, inclusive.
+
+GraphBase programs are designed to produce identical results on almost
+all existing computers and operating systems. An improved version of the
+portable subtractive method recommended in {\sl Seminumerical Algorithms},
+Section~3.6, is used to generate random numbers in the routines below.
+The period length of the generated numbers is at least $2^{55}-1$, and
+it is in fact plausibly conjectured to be $2^{85}-2^{30}$ for all but
+at most one choice of the |seed| value. The low-order bits of the
+generated numbers are just as random as the high-order bits.
+
+@ Changes might be needed when these routines are ported to different
+systems, because the programs have been written to be most efficient
+on binary computers that use two's complement notation. Almost all
+modern computers are based on two's complement arithmetic, but if you have a
+nonconformist machine you may have to revise the code in sections that
+are listed under `system dependencies' in the index.
+
+A validation program is provided so that installers can tell if |gb_flip|
+is working properly. To make the test, simply run |test_flip|.
+
+@(test_flip.c@>=
+#include <stdio.h>
+#include "gb_flip.h" /* all users of |gb_flip| should do this */
+main()
+{@+int j;
+ gb_init_rand(-314159);
+ if (gb_next_rand()!=119318998) {
+ fprintf(stderr,"Failure on the first try!\n"); return -1;
+ }
+ for (j=1; j<=133; j++)
+ gb_next_rand();
+ if (gb_unif_rand(0x55555555)!=748103812) {
+ fprintf(stderr,"Failure on the second try!\n"); return -2;
+ }
+ fprintf(stderr,"OK, the gb_flip routines seem to work!\n");
+}
+
+@ The \Cee\ code for |gb_flip| doesn't have a main routine; it's just a
+bunch of subroutines to be incorporated into programs at a higher level,
+via the system loading routine. Here is the general outline of \.{gb\_flip.c}:
+
+@p
+@<Private declarations@>@;
+@<External declarations@>@;
+@<External functions@>
+
+@* The subtractive method. If $m$ is any even number, and if the
+numbers $a_0$, $a_1$, \dots,~$a_{54}$ are not all even, then the numbers
+generated by the recurrence
+$$ a_n=(a_{n-24}-a_{n-55})\bmod m $$
+have a period length of at least $2^{55}-1$, because the residues
+$a_n\bmod2$ have a period of this length. Furthermore, the numbers 24 and~55
+in this recurrence are sufficiently large that deficiencies in randomness
+due to the simplicity of the recurrence are negligible in most applications.
+
+Here we take $m=2^{31}$ so as to get the full set of nonnegative numbers
+on a 32-bit computer. The recurrence is computed by maintaining an array
+of 55 values, $A[1]\ldots A[55]$. We also set |A[0]=-1| to act as a sentinel.
+
+@<Private...@>=
+static long A[56] = {-1}; /* pseudo-random values */
+
+@ Every external variable should be declared twice in this \.{CWEB} file:
+once for |gb_flip| itself (the ``real'' declaration for storage allocation
+purposes), and once in \.{gb\_flip.h} (for cross-references by |gb_flip|
+users).
+
+The pointer variable |gb_flip_ptr| should not actually be mentioned explicitly
+by user routines; it is made public only for efficiency, so that the
+|gb_next_rand| macro can access the private |A| table.
+
+@<External declarations@>=
+long *gb_flip_ptr=A; /* the next |A| value to be exported */
+
+@ Incidentally, we hope that optimizing compilers are smart enough to
+do the right thing with |gb_next_rand|.
+
+@d gb_next_rand() (*gb_flip_ptr>=0? *gb_flip_ptr--: gb_flip_cycle())
+
+@(gb_flip.h@>=
+#define gb_next_rand()@t\quad@>(*gb_flip_ptr>=0?*gb_flip_ptr--:gb_flip_cycle())
+extern long *gb_flip_ptr; /* the next |A| value to be used */
+extern long gb_flip_cycle(); /* compute 55 more pseudo-random numbers */
+
+@ The user is not supposed to call |gb_flip_cycle| directly either.
+It is a routine invoked by the macro |gb_next_rand()| when |gb_flip_ptr|
+points to the negative value in |A[0]|.
+
+The purpose of |gb_flip_cycle| is to do 55 more steps of the basic
+recurrence, at high speed, and to reset |gb_flip_ptr|.
+
+The nonnegative remainder of $(x-y)\bmod 2^{31}$ is computed here by
+doing a logical-and with the constant |0x7fffffff|. On computers without
+two's complement arithmetic it may be more efficient to add the
+value $2^{30}$ twice to $(x-y)$, if $(x-y)$ turns out to be negative.
+Careful calculations are essential to preserve system-independence of
+the GraphBase results.
+@^system dependencies@>
+
+The sequence of random numbers returned by successive calls of |gb_next_rand()|
+isn't really $a_n$, $a_{n+1}$, \dots, as defined by the basic recurrence above;
+blocks of 55 consecutive values are essentially being ``flipped'' or
+``reflected,'' i.e., output in reverse order, because |gb_next_rand()|
+makes the value of |gb_flip_ptr| decrease instead of increase.
+But such flips don't make the results any less random.
+
+@<External functions@>=
+long gb_flip_cycle()
+{@+register long *ii, *jj;
+ for (ii=&A[1],jj=&A[32];jj<=&A[55];ii++,jj++)
+ *ii=(*ii-*jj)&0x7fffffff;
+ for (jj=&A[1];ii<=&A[55];ii++,jj++)
+ *ii=(*ii-*jj)&0x7fffffff;
+ gb_flip_ptr=&A[54];
+ return A[55];
+}
+
+@* Initialization. To get everything going, we use a scheme like that
+recommended in {\sl Seminumerical Algorithms}, but revised so that the
+least significant bits of the starting values depend on the entire
+seed, not just on the seed's least significant bits.
+
+Notice that we jump around in the array by increments of 21, a number that is
+relatively prime to~55. Repeated skipping by steps of 21~mod~55 keeps the
+values we're computing spread out as far from each other as possible in the
+array, since 21, 34, and 55 are consecutive
+Fibonacci numbers (see the discussion of Fibonacci hashing in
+Section 6.4 of {\sl Sorting and Searching\/}). Our initialization mechanism
+would be rather poor if we didn't do something like that to disperse the values
+(see {\sl Seminumerical Algorithms}, exercise 3.2.2--2).
+
+@<External f...@>=
+void gb_init_rand(seed)
+ long seed;
+{@+register int i;
+ register long prev=seed, next=1;
+ seed=prev=prev & 0x7fffffff; /* strip off the sign */
+ A[55]=prev;
+ for (i=21; i; i=(i+21)%55) {
+ A[i]=next;
+ @<Compute a new |next| value, based on |next|, |prev|, and |seed|@>;
+ prev=A[i];
+ }
+ @<Get the array values ``warmed up''@>;
+}
+
+@ Here we have two more instances of $(x-y)\bmod 2^{31}$ that should
+be computed in some other way on atypical machines.
+@^system dependencies@>
+
+Incidentally, if |test_flip| fails, the person debugging these routines will
+want to know some of the intermediate numbers computed during initialization.
+The first nontrivial values calculated by |gb_init_rand| are
+|A[42]=2147326568|; |A[8]=1073977445|; |A[29]=536517481|.
+Once you get these right, the rest should be easy.
+
+An early version of this routine simply said `|seed>>1|' instead of making
+|seed| shift cyclically. This method had an interesting flaw:
+When the original |seed| was a number of the form $4s+1$, the first
+54 elements $A[1]$, \dots,~$A[54]$ were set to exactly the same values
+as when |seed| was $4s+2$. Therefore one out of every four seed values
+was effectively being wasted.
+
+@<Compute a new |next|...@>=
+next=(prev-next) & 0x7fffffff;
+if (seed&1) seed=0x40000000+(seed>>1);
+else seed>>=1; /* cyclic shift right 1 */
+next=(next-seed) & 0x7fffffff;
+
+@ After the first 55 values have been computed as a function of |seed|,
+they aren't random enough for us to start using them right away. For example,
+we have set |A[21]=1|, in order to ensure that at least one starting value
+is an odd number. But once the sequence $a_n$ gets going far enough from
+its roots, the initial transients become imperceptible. Therefore we will call
+|gb_flip_cycle| five times, effectively skipping past the first 275
+elements of the sequence; this has the desired effect, and it also
+initializes |gb_flip_ptr|.
+
+Note: It is possible to express the least significant bit of the
+generated numbers as a linear combination mod~2 of the 31 bits of
+|seed| and of the constant~1. For example, the first generated number
+turns out to be odd if and only if
+$$s_{24}+s_{23}+s_{22}+s_{21}+s_{19}+s_{18}+s_{15}+s_{14}+s_{13}+s_{11}+
+s_{10}+s_{8}+s_{7}+s_{6}+s_{2}+s_{1}+s_{0}$$ is odd, when
+$|seed|=(s_{31}\ldots s_1s_0)_2$. We can represent this linear
+combination conveniently by the hexadecimal number |0x01ecedc7|; the
+\.1 stands for $s_{24}$ and the final \.7 stands for $s_2+s_1+s_0$.
+The first ten least-significant bits turn out to be respectively
+|0x01ecedc7|, |0xdbbdc362|, |0x400e0b06|, |0x0eb73780|, |0xda0d66ae|,
+|0x002b63bc|, |0xadb801ed|, |0x8077bbbc|, |0x803d9db5|, and
+|0x401a0eda| in this notation (using the sign bit to indicate cases
+when 1 must be added to the sum).
+
+We must admit that these ten 32-bit patterns do not look at all
+random; the number of \.b's, \.d's, and \.0's is unusually high. (Before
+the ``warmup cycles,'' the patterns are even more regular.) This
+phenomenon eventually disappears, however, as the sequence proceeds;
+and it does not seem to imply any serious deficiency in practice, even
+at the beginning of the sequence, once we've done the warmup exercises.
+
+@<Get the array...@>=
+(void) gb_flip_cycle();
+(void) gb_flip_cycle();
+(void) gb_flip_cycle();
+(void) gb_flip_cycle();
+(void) gb_flip_cycle();
+
+@ @(gb_flip.h@>=
+extern void gb_init_rand();
+
+@* Uniform integers.
+Here is a simple routine that produces a uniform integer between
+0 and~$m-1$, inclusive, when $m$ is any positive integer less than $2^{31}$.
+It avoids the bias toward small values that would occur if we simply
+calculated |gb_next_rand()%m|. (The bias is insignificant when |m| is
+small, but it can be serious when |m| is large. For example, if
+$m\approx 2^{32}/3$, the simple remainder algorithm would give an answer
+less than $m/2$ about 2/3 of the time.)
+
+This routine consumes fewer than two random numbers, on the average,
+for any fixed~$m$.
+
+In the |test_flip| program this routine should compute |t=m|,
+then it should reject the values |r=2081307921|, 1621414801, and
+1469108743 before returning the answer 748103812.
+
+@d two_to_the_31 ((unsigned long)0x80000000)
+
+@<External f...@>=
+long gb_unif_rand(m)
+ long m;
+{@+register unsigned long t=two_to_the_31-(two_to_the_31 % m);
+ register long r;
+ do {
+ r=gb_next_rand();
+ } while (t<=(unsigned long)r);
+ return r%m;
+}
+
+@ @(gb_flip.h@>=
+extern long gb_unif_rand();
+
+@* Index. Here is a list that shows where the identifiers of this program are
+defined and used.