\input cwebmac
% This file is part of the Stanford GraphBase (c) Stanford University 1993
% This material goes at the beginning of all Stanford GraphBase CWEB files

\def\topofcontents{
  \leftline{\sc\today\ at \hours}\bigskip\bigskip
  \centerline{\titlefont\title}}

\font\ninett=cmtt9
\def\botofcontents{\vskip 0pt plus 1filll
    \ninerm\baselineskip10pt
    \noindent\copyright\ 1993 Stanford University
    \bigskip\noindent
    This file may be freely copied and distributed, provided that
    no changes whatsoever are made. All users are asked to help keep
    the Stanford GraphBase files consistent and ``uncorrupted,''
    identical everywhere in the world. Changes are permissible only
    if the modified file is given a new name, different from the names of
    existing files in the Stanford GraphBase, and only if the modified file is
    clearly identified as not being part of that GraphBase.
    (The {\ninett CWEB} system has a ``change file'' facility by
    which users can easily make minor alterations without modifying
    the master source files in any way. Everybody is supposed to use
    change files instead of changing the files.)
    The author has tried his best to produce correct and useful programs,
    in order to help promote computer science research,
    but no warranty of any kind should be assumed.
    \smallskip\noindent
    Preliminary work on the Stanford GraphBase project
    was supported in part by National Science
    Foundation grant CCR-86-10181.}

\def\prerequisite#1{\def\startsection{\noindent
    Important: Before reading {\sc\title},
    please read or at least skim the program for {\sc#1}.\bigskip
    \let\startsection=\stsec\stsec}}
\def\prerequisites#1#2{\def\startsection{\noindent
    Important: Before reading {\sc\title}, please read
    or at least skim the programs for {\sc#1} and {\sc#2}.\bigskip
    \let\startsection=\stsec\stsec}}

\def\title{GB\_\,FLIP}
\hyphenation{semi-numer-ical}


\N{1}{1}Introduction. This is {\sc GB\_\,FLIP}, the module used by GraphBase
programs to generate random numbers.

To use the routines in this file, first call the function \PB{\\{gb\_init%
\_rand}(\\{seed})}.
Subsequent uses of the macro \PB{\\{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 \PB{\\{seed}} value.  The low-order bits of the
generated numbers are just as random as the high-order bits.

\fi

\M{2}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 might 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
{\sc GB\_\,FLIP} is working properly. To make the test, simply run
\.{test\_flip}.

\Y\B\4\X2:\.{test\_flip.c }\X${}\E{}$\6
\8\#\&{include} \.{<stdio.h>}\6
\8\#\&{include} \.{"gb\_flip.h"}\C{ all users of {\sc GB\_\,FLIP} should do
this }\7
\&{int} \\{main}(\,)\1\1\2\2\6
${}\{{}$\5
\1\&{long} \|j;\7
${}\\{gb\_init\_rand}({-}\T{314159\$L});{}$\6
\&{if} ${}(\\{gb\_next\_rand}(\,)\I\T{119318998}){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"Failure\ on\ the\ firs}\)\.{t\ try!\\n"});{}$\6
\&{return} ${}{-}\T{1};{}$\6
\4${}\}{}$\2\6
\&{for} ${}(\|j\K\T{1};{}$ ${}\|j\Z\T{133};{}$ ${}\|j\PP){}$\1\5
\\{gb\_next\_rand}(\,);\2\6
\&{if} ${}(\\{gb\_unif\_rand}(\T{\^55555555\$L})\I\T{748103812}){}$\5
${}\{{}$\1\6
${}\\{fprintf}(\\{stderr},\39\.{"Failure\ on\ the\ seco}\)\.{nd\ try!\\n"});{}$%
\6
\&{return} ${}{-}\T{2};{}$\6
\4${}\}{}$\2\6
${}\\{fprintf}(\\{stderr},\39\.{"OK,\ the\ gb\_flip\ rou}\)\.{tines\ seem\ to\
work!\\}\)\.{n"});{}$\6
\&{return} \T{0};\6
\4${}\}{}$\2\par
\fi

\M{3}The \CEE/ code for {\sc 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}:

\Y\B\X4:Private declarations\X\6
\X5:External declarations\X\6
\X7:External functions\X\par
\fi

\N{1}{4}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-55}-a_{n-24})\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 that we 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 \PB{$\|A[\T{0}]\K{-}\T{1}$} to
act as a sentinel.

\Y\B\4\X4:Private declarations\X${}\E{}$\6
\&{static} \&{long} \|A[\T{56}]${}\K\{{-}\T{1}\}{}$;\C{ pseudo-random values }%
\par
\U3.\fi

\M{5}Every external variable should be declared twice in this \.{CWEB}
file: once for {\sc GB\_\,FLIP} itself (the ``real'' declaration for
storage allocation purposes), and once in \.{gb\_flip.h} (for
cross-references by {\sc GB\_\,FLIP} users).

The pointer variable \PB{\\{gb\_fptr}} should not be mentioned explicitly
by user routines. It is made public only for efficiency, so that the
\PB{\\{gb\_next\_rand}} macro can access the private \PB{\|A} table.

\Y\B\4\X5:External declarations\X${}\E{}$\6
\&{long} ${}{*}\\{gb\_fptr}\K\|A{}$;\C{ the next \PB{\|A} value to be exported
}\par
\U3.\fi

\M{6}The numbers generated by \PB{\\{gb\_next\_rand}(\,)} seem to be
satisfactory for most
purposes, but they do fail a stringent test called the ``birthday spacings
test,'' devised by George Marsaglia. [See, for example, {\sl Statistics and
Probability Letters\/ \bf8} (1990), 35--39.] One way to get numbers that
pass the birthday test is to discard half of the values, for example
by changing `\PB{\\{gb\_flip\_cycle}(\,)}' to `\PB{$(\\{gb\_flip\_cycle}(\,),%
\\{gb\_flip\_cycle}(\,))$}'
in the definition of \PB{\\{gb\_next\_rand}(\,)}. Users who wish to make such a
change
should define their own substitute macro.

Incidentally, we hope that optimizing compilers are smart enough to
do the right thing with \PB{\\{gb\_next\_rand}}.

\Y\B\4\D$\\{gb\_next\_rand}()$ \5
$({*}\\{gb\_fptr}\G\T{0}\?{*}\\{gb\_fptr}\MM:\\{gb\_flip\_cycle}(\,){}$)\par
\Y\B\4\X6:\.{gb\_flip.h }\X${}\E{}$\6
\8\#\&{define} \\{gb\_next\_rand}(\,) \5${}\hbox{\quad}({*}\\{gb\_fptr}\G\T{0}%
\?{*}\\{gb\_fptr}\MM:\\{gb\_flip\_cycle}(\,)){}$\6
\&{extern} \&{long} ${}{*}\\{gb\_fptr}{}$;\C{ the next \PB{\|A} value to be
used }\6
\&{extern} \&{long} \\{gb\_flip\_cycle}(\,);\C{ compute 55 more pseudo-random
numbers }\par
\As11\ET13.\fi

\M{7}The user is not supposed to call \PB{\\{gb\_flip\_cycle}} directly either.
It is a routine invoked by the macro \PB{\\{gb\_next\_rand}(\,)} when \PB{\\{gb%
\_fptr}}
points to the negative value in \PB{\|A[\T{0}]}.

The purpose of \PB{\\{gb\_flip\_cycle}} is to do 55 more steps of the basic
recurrence, at high speed, and to reset \PB{\\{gb\_fptr}}.

The nonnegative remainder of $(x-y)\bmod 2^{31}$ is computed here by
doing a logical-and with the constant \PB{\T{\^7fffffff}}. This technique
doesn't work on computers that do not perform two's complement
arithmetic. An alternative for such machines is to add the value
$2^{30}$ twice to $(x-y)$, when $(x-y)$ turns out to be negative.
Careful calculations are essential because the GraphBase results
must be identical on all computer systems.

The sequence of random numbers returned by successive calls of \PB{\\{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''---output in reverse order---because \PB{\\{gb\_next\_rand}(\,)}
makes the value of \PB{\\{gb\_fptr}} decrease instead of increase.
But such flips don't make the results any less random.

\Y\B\4\D$\\{mod\_diff}(\|x,\|y)$ \5
$(((\|x)-(\|y))\AND\T{\^7fffffff}{}$)\C{ difference modulo $2^{31}$ }\par
\Y\B\4\X7:External functions\X${}\E{}$\6
\&{long} \\{gb\_flip\_cycle}(\,)\1\1\2\2\6
${}\{{}$\5
\1\&{register} \&{long} ${}{*}\\{ii},{}$ ${}{*}\\{jj};{}$\7
\&{for} ${}(\\{ii}\K{\AND}\|A[\T{1}],\39\\{jj}\K{\AND}\|A[\T{32}];{}$ ${}\\{jj}%
\Z{\AND}\|A[\T{55}];{}$ ${}\\{ii}\PP,\39\\{jj}\PP){}$\1\5
${}{*}\\{ii}\K\\{mod\_diff}({*}\\{ii},\39{*}\\{jj});{}$\2\6
\&{for} ${}(\\{jj}\K{\AND}\|A[\T{1}];{}$ ${}\\{ii}\Z{\AND}\|A[\T{55}];{}$ ${}%
\\{ii}\PP,\39\\{jj}\PP){}$\1\5
${}{*}\\{ii}\K\\{mod\_diff}({*}\\{ii},\39{*}\\{jj});{}$\2\6
${}\\{gb\_fptr}\K{\AND}\|A[\T{54}];{}$\6
\&{return} \|A[\T{55}];\6
\4${}\}{}$\2\par
\As8\ET12.
\U3.\fi

\N{1}{8}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).

\Y\B\4\X7:External functions\X${}\mathrel+\E{}$\6
\&{void} \\{gb\_init\_rand}(\\{seed})\1\1\6
\&{long} \\{seed};\2\2\6
${}\{{}$\5
\1\&{register} \&{long} \|i;\6
\&{register} \&{long} \\{prev}${}\K\\{seed},{}$ \\{next}${}\K\T{1};{}$\7
${}\\{seed}\K\\{prev}\K\\{mod\_diff}(\\{prev},\39\T{0}){}$;\C{ strip off the
sign }\6
${}\|A[\T{55}]\K\\{prev};{}$\6
\&{for} ${}(\|i\K\T{21};{}$ \|i; ${}\|i\K(\|i+\T{21})\MOD\T{55}){}$\5
${}\{{}$\1\6
${}\|A[\|i]\K\\{next};{}$\6
\X9:Compute a new \PB{\\{next}} value, based on \PB{\\{next}}, \PB{\\{prev}},
and \PB{\\{seed}}\X;\6
${}\\{prev}\K\|A[\|i];{}$\6
\4${}\}{}$\2\6
\X10:Get the array values ``warmed up''\X;\6
\4${}\}{}$\2\par
\fi

\M{9}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
\PB{\\{gb\_init\_rand}} are \PB{$\|A[\T{42}]\K\T{2147326568}$}, \PB{$\|A[\T{8}]%
\K\T{1073977445}$}, and
\PB{$\|A[\T{29}]\K\T{536517481}$}.  Once you get those right, the rest should
be easy.

An early version of this routine simply said `\PB{$\\{seed}\GG\T{1}$}' instead
of making
\PB{\\{seed}} shift cyclically. This method had an interesting flaw:
When the original \PB{\\{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 \PB{\\{seed}} was $4s+2$. Therefore one out of every four seed values
was effectively being wasted.

\Y\B\4\X9:Compute a new \PB{\\{next}} value, based on \PB{\\{next}}, \PB{%
\\{prev}}, and \PB{\\{seed}}\X${}\E{}$\6
$\\{next}\K\\{mod\_diff}(\\{prev},\39\\{next});{}$\6
\&{if} ${}(\\{seed}\AND\T{1}){}$\1\5
${}\\{seed}\K\T{\^40000000}+(\\{seed}\GG\T{1});{}$\2\6
\&{else}\1\5
${}\\{seed}\MRL{{\GG}{\K}}\T{1}{}$;\C{ cyclic shift right 1 }\2\6
${}\\{next}\K\\{mod\_diff}(\\{next},\39\\{seed}){}$;\par
\U8.\fi

\M{10}After the first 55 values have been computed as a function of \PB{%
\\{seed}},
they aren't random enough for us to start using them right away. For example,
we have set \PB{$\|A[\T{21}]\K\T{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 call
\PB{\\{gb\_flip\_cycle}} five times, effectively skipping past the first 275
elements of the sequence; this has the desired effect. It also
initializes \PB{\\{gb\_fptr}}.

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
\PB{\\{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
$\PB{\\{seed}}=(s_{31}\ldots s_1s_0)_2$.  We can represent this linear
combination conveniently by the hexadecimal number \PB{\T{\^01ecedc7}}; 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
\PB{\T{\^01ecedc7}}, \PB{\T{\^dbbdc362}}, \PB{\T{\^400e0b06}}, \PB{\T{%
\^0eb73780}}, \PB{\T{\^da0d66ae}},
\PB{\T{\^002b63bc}}, \PB{\T{\^adb801ed}}, \PB{\T{\^8077bbbc}}, \PB{\T{%
\^803d9db5}}, and
\PB{\T{\^401a0eda}} 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.

\Y\B\4\X10:Get the array values ``warmed up''\X${}\E{}$\6
(\&{void}) \\{gb\_flip\_cycle}(\,);\6
(\&{void}) \\{gb\_flip\_cycle}(\,);\6
(\&{void}) \\{gb\_flip\_cycle}(\,);\6
(\&{void}) \\{gb\_flip\_cycle}(\,);\6
(\&{void}) \\{gb\_flip\_cycle}(\,);\par
\U8.\fi

\M{11}\B\X6:\.{gb\_flip.h }\X${}\mathrel+\E{}$\6
\&{extern} \&{void} \\{gb\_init\_rand}(\,);\par
\fi

\N{1}{12}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 \PB{$\\{gb\_next\_rand}(\,)\MOD\|m$}. (The bias is insignificant
when \PB{\|m} is
small, but it can be serious when \PB{\|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 (\PB{\\{main}}), this routine should compute %
\PB{$\|t\K\|m$},
then it should reject the values \PB{$\|r\K\T{2081307921}$}, 1621414801, and
1469108743 before returning the answer 748103812.

\Y\B\4\D$\\{two\_to\_the\_31}$ \5
((\&{unsigned} \&{long}) \T{\^80000000})\par
\Y\B\4\X7:External functions\X${}\mathrel+\E{}$\6
\&{long} \\{gb\_unif\_rand}(\|m)\1\1\6
\&{long} \|m;\2\2\6
${}\{{}$\5
\1\&{register} \&{unsigned} \&{long} \|t${}\K\\{two\_to\_the\_31}-(\\{two\_to%
\_the\_31}\MOD\|m);{}$\6
\&{register} \&{long} \|r;\7
\&{do}\5
${}\{{}$\1\6
${}\|r\K\\{gb\_next\_rand}(\,);{}$\6
\4${}\}{}$\5
\2\5
\&{while} ${}(\|t\Z{}$(\&{unsigned} \&{long}) \|r);\6
\&{return} ${}\|r\MOD\|m;{}$\6
\4${}\}{}$\2\par
\fi

\M{13}\B\X6:\.{gb\_flip.h }\X${}\mathrel+\E{}$\6
\&{extern} \&{long} \\{gb\_unif\_rand}(\,);\par
\fi

\N{1}{14}Index. Here is a list that shows where the identifiers of this program
are
defined and used.
\fi


\inx
\fin
\con
