/*
Copyright 1995 David M. Dantowitz
All Rights Reserved.
LEGAL STUFF:
Non-commercial distribution and non-commercial use of this code
are permitted, provided the file is distributed in its original
form. Minor modifications are permitted but this notice and the
documentation included with the code may not be removed.
The decision to use the code, methods, and algorithms expressed
below is done at the discresion of the user. The author has
made a serious attempt to verify that this code works as
described and has provided substantial diagnostic tests, but
makes no claims or warranties to that effect.
If distributed or used commercially, modified, or used in any form
that requires substantial modifications (e.g., translation or
re-coding in microcode, assembly language, etc.), please contact
the author so he may coordinate any such efforts and provide proper
licensing information.
WHAT THE CODE DOES:
This code generates random numbers using very fast MOD 2^N
arithmetic. The implementation below uses only generic C
code and nothing processor or compiler dependent. The code
offers a method to compute X(n+1) from Xn and X(n+i) from
Xn for any non-negative i (e.g., skipping ahead i numbers
in a sequence).
See the routine SampleUse() for an example of how to use this
generator. The runTest() routine will test to see if the code
is properly configured to work with your compiler and machine
architecture. It is VERY important that you run this routine
to verify that the math is being performed properly. If not
you could be using VERY UNRANDOM numbers for your work.
Due to the nature of the computation, substantial diagnostics
are provided to ensure that the computations are performed
correctly. This project found 3 compiler bugs in two compilers
on two different platforms. You may have a similar experience.
If the runTest diagnostics find a problem, look for the keyword
"HINT" in the code for areas of the code you can re-configure
to run properly on your system.
This code is an out-growth of my research on parallel simulation.
With the need to generate random numbers in a parallel environment
I came across the NAS Benchmarks' method of computing random
numbers. Since the generator used in the NAS Benchmarks uses
mod 2^46 arithmetic, I set out with the goal of computing mod 2^N
arithmetic as fast as possible using generic C code. Below you
will find my results on this subject.
If you have received the mini version of this code, and would like
the larger tutorial version, please send the author a note.
RESUME/ADVERTISEMENT:
If you find that using this code improves your performance, please
drop me a line to let me know. If you need any consulting expertise
in this area I'd love to hear from you. I'm currently completing
my M.S. thesis research and will be looking for some consulting
positions in the near future. Also, if you have any spare cycles
on a shared memory machine, I'm always looking to port my simulations
to new machines and run some benchmarks.
I also have a substantial background in processor verification and
automated testing (this self-testing code is a simplified version of
that type of work). If you need to find processor bugs (from FP
testing to pipeline logic) or need to automate a testing process
let me know. (I've tested everything from software systems,
compilers, automated teller machines, and peripherals, to new CISC
and RISC machines.).
And just for fun, my other area of interest is multimedia CD-ROMs.
Let me know if you need assistance creating CD-ROMs (multimedia,
audio CDs or otherwise).
If you have suggestions, or need assistance in the use, porting,
licensing, reporting bugs, or questions on optimizing other
computations or software development, feel free to contact me at:
Dantowitz Consulting & Research
dantowitz@dantowitz.com
David Dantowitz
---------------------------------------------------------------------
BACKGROUND:
Most of the code in this file serves as a test bed and tutorial on
performing fast modulo 2^N arithmetic with generic C code.
The core motivation is to compute the series: X , X , X , X , ...
0 1 2 3
where:
N
X = AX mod 2
k+1 k
This series forms the basis of a random number generator. See the
routine SampleUse() for an example of how to use this generator.
Code exists to compute a sequence of values or to skip ahead ANY
number of terms in a sequence (you may even skip ahead in increments
greater than 2^32!)
Besides speed, the major goal here was to use generic C code and
avoid machine specific code. If faster methods are created using
assembly language, the generic versions in C may serve as a source
for testing.
The random number generator used here has a period on the order of
2^44 (1.76x10^13). It uses a linear congruential recursion with
the values A=5^13 and N=46 (or 52). This generator is based on
information from "The NAS Parallel Benchmarks," by D. Bailey, et al,
RNR Technical Report RNR-94-007, March 1994. For more information
see http://www.nas.nasa.gov.
The same paper says:
"the radix, 2^46 may be changed to 2^52 if all systems properly
implement IEEE-754. If this is done, the generator's period
increases to approximately 1.13x10^15."
To make this change, all you have to do is change USE_2_46 from
1 to 0. This will change the radix from 2^46 to 2^52.
HINT: This code takes advantage of the representation of doubles to
manipulate the exponent and mantissa directly. If doubles are not
8 bytes long on your system, send me mail and I will send you code
to work with 10-byte or 12-byte doubles. On the systems I tested,
using 8-byte doubles was preferable to using 10- or 12-byte doubles.
Contact David Dantowitz at dantowitz@aol.com for help with anything.
David Dantowitz April 18, 1995 Version 1.7
. Initial external distribution.
David Dantowitz March 23, 1995 Version 1.6
. Fix code to work around a bug in another compiler.
David Dantowitz March 20, 1995 Version 1.5
. Final touches for first external release.
David Dantowitz March 14, 1995 Version 1.4
. Find bug in one optimizer and add #define BUGGY_OPTIMIZER to
handle the problem.
David Dantowitz Febrary 28, 1995 Version 1.3
. Improve performance.
David Dantowitz Febrary 21, 1995 Version 1.2
. Clean up and add lots of self-testing code.
. Add a routine to skip n calls to the generator.
. Add code to work with 8-byte IEEE doubles.
. Removed 10 & 12 byte long double code.
David Dantowitz January 21, 1995 Version 1.1
. Add some helpful documentation.
David Dantowitz January 7, 1995 Version 1.0
. Initial implementation.
*/
#include
#include
#include
#include
/*
This controls one small bit of code in the skipNRandCalls0
routine. At least one compiler had a bug in its optimizer.
Set to 1 if you get the error:
"the skip N seeds routine fails for very large N"
If this does not fix the problem, please contact the author.
*/
#define BUGGY_OPTIMIZER 0
/*
Set STANDALONE to 1 in order to run the test and performance
software.
Set STANDALONE to 0 in order to use these routines from other
modules.
*/
#define STANDALONE 1
/*
The diagnostic tests will verify that this is set properly.
*/
#define BIGENDIAN 1
/*
Choose routines or macros (in-line expansion)
(The routine versions are useful when debugging new compilers.)
*/
#define PARTSROUTINE 0
#if PARTSROUTINE
#define parts PartsR
#else
#define parts PartsM
#endif
#if PARTSROUTINE
#define parts2 PartsR2
#else
#define parts2 PartsM2
#endif
/*
Use 5^13 as the multiplier.
*/
#define MULTIPLIER 1220703125
#define INITIAL_SEED 314159265
/*
HINT: Please ensure that the EightByteDouble type uses 8 bytes
of storage. (This, and many other features are tested when compiled
with STANDALONE set to 1.)
*/
typedef double EightByteDouble;
/*
This typedef must be a four byte integer.
*/
typedef signed long int FourByteInt;
typedef unsigned long int UFourByteInt;
/*
Set to 1 for mod 2^46 arithmetic.
Set to 0 for mod 2^52 arithmetic.
*/
#define USE_2_46 1
#if USE_2_46
#define FIXED_POWER 46
/*
These are the values for the 100K-th seed for an initial
seed of 314159265. Note that the mantissae are stored
using a 64-bit representation.
*/
#define SEED_100K_LO 0xC2840000
#define SEED_100K_HI 0x95C5FAA1
#define SEED_100K_M 45
#else
#define FIXED_POWER 52
/*
These are the values for the 100K-th seed for an initial
seed of 314159265. Note that the mantissae are stored
using a 64-bit representation.
*/
#define SEED_100K_LO 0x870A1000
#define SEED_100K_HI 0x8E5717EA
#define SEED_100K_M 51
#endif
/*
This is one half of FIXED_POWER.
*/
#define FIXED_POWER_HALF (FIXED_POWER >> 1)
/*
Divide the multiplier by 2^FIXED_POWER_HALF
e.g., for FIXED_POWER_HALF = 23, we have:
A = 5^13 = 1,220,703,125 = 145 * 2^23 + 4,354,965
^ ^
MULT_HI MULT_LO
*/
#define MULT_HI (MULTIPLIER >> FIXED_POWER_HALF)
#define MULT_LO (MULTIPLIER & ( (1L<> 21)
/*
KEEP_N_PLUS_1_BITS will be asked to keep from 1 to 26 bits from
a mantissa with 1 to 52 non-zero bits. The high N bits will be
kept and the lower bits will be masked out.
In other words, this code is performing the following operation:
FPout = (Floor(FPin/2^P)) * 2^P
HINT:
If you have problems running the program, change to a different
version of KEEP_N_PLUS_1_BITS. Some compilers have problems
shifting values properly, others may be faster using one
variant of this core macro than another.
At least one compiler has a bug shifting more than 47 bits.
As a result, you should use a version of KEEP_N_PLUS_1_BITS that
never tries to shift more than 47 bits. (See the one that uses
(27 - N) below.) Others may have optimizer bugs.
You may also have to change: USE_SMALL_VERSION if you still have
bugs.
*/
#if 1
#define KEEP_N_PLUS_1_BITS(FPin, FPout, N) \
{\
(FPout).u.fpParts.lo = \
(FPin).u.fpParts.lo & (0xFFFFFFFF << (52 - (N))); \
\
(FPout).u.fpParts.hi = \
(FPin).u.fpParts.hi & (((FourByteInt) 0xFFF00000) >> (N)); \
}
#endif
#if 0
#define KEEP_N_PLUS_1_BITS(FPin, FPout, N) \
{\
(FPout).u.fpParts.lo = \
(FPin).u.fpParts.lo & (0xFE000000 << (27 - (N))); \
\
(FPout).u.fpParts.hi = \
(FPin).u.fpParts.hi & (((FourByteInt) 0xFFF00000) >> (N)); \
}
#endif
#if 0
#define KEEP_N_PLUS_1_BITS(FPin, FPout, N) \
{\
(FPout).u.fpParts.lo = \
(FPin).u.fpParts.lo & (0xFFFFFFFF << (52 - (N))); \
\
(FPout).u.fpParts.hi = \
N<20 ? \
((FPin).u.fpParts.hi & (0xFFFFFFFF << (20 - (N)))) : \
(FPin).u.fpParts.hi; \
}
#endif
#if 0
#define KEEP_N_PLUS_1_BITS(FPin, FPout, N) \
{\
register int * MASK = &(mantissaMasks[N << 1]);\
\
(FPout).u.fpParts.lo = (FPin).u.fpParts.lo & MASK[0]; \
\
(FPout).u.fpParts.hi = (FPin).u.fpParts.hi & MASK[1]; \
}
#endif
#if 0
#define KEEP_N_PLUS_1_BITS(FPin, FPout, N) \
KEEP_N_PLUS_1_BITSh(&FPin, &FPout, N);
void KEEP_N_PLUS_1_BITSh(fp *FPin, fp *FPout, int N);
void
KEEP_N_PLUS_1_BITSh(fp *FPin, fp *FPout, int N)
{
FPout->u.fpParts.lo =
FPin->u.fpParts.lo & (0xFE000000 << (27 - (N)));
/*
FPout->u.fpParts.lo =
FPin->u.fpParts.lo & (0xFFFFFFFF << (52 - (N)));
*/
FPout->u.fpParts.hi =
FPin->u.fpParts.hi & (((FourByteInt) 0xFFF00000) >> (N));
}
#endif
/*
HINT:
Set to 1 for a slight optimization on most machines. If you
experience any problems, you may want to set it to 0.
*/
#define USE_SMALL_VERSION 1
#if USE_SMALL_VERSION
/*
This version is used when we know we will never have N greater than
20. (See mod2Nf, which, intern, is used by nextRandFinal.)
*/
#define KEEP_N_PLUS_1_BITS_SMALL(FPin, FPout, N) \
{\
(FPout).u.fpParts.lo = 0; \
\
(FPout).u.fpParts.hi = \
(FPin).u.fpParts.hi & (((FourByteInt) 0xFFF00000) >> (N)); \
}
#else
#define KEEP_N_PLUS_1_BITS_SMALL(FPin, FPout, N) \
KEEP_N_PLUS_1_BITS(FPin,FPout,N)
#endif
/*
Table of contents (routines).
*/
#if PARTSROUTINE
/*
Utility routines (replaced by macros for performance.)
*/
void PartsR(fp *dividend_residual, fp *quotient, int power);
void PartsR2(fp *dividend_residual, fp *quotient, int power);
#endif
/*
Generic mod P multiplication
(computes c = (a * b) mod 2^FIXED_POWER).
*/
void mulABmodP(fp a, fp b, fp *c);
/*
Computes c = (a * a) mod 2^FIXED_POWER.
*/
void squareModP(fp a, fp *c);
/*
Set the random seed.
*/
void setRandSeed(EightByteDouble s);
/*
Get the current random seed.
*/
EightByteDouble getRandSeed(void);
/*
Three fast random number routines.
*/
EightByteDouble nextRandEarly(void);
EightByteDouble nextRandSecond(void);
EightByteDouble nextRandFinal(void);
/*
A routine to skip ahead in a sequence.
*/
EightByteDouble
skipNRandCalls0(EightByteDouble initialSeed, EightByteDouble nFP);
/*
Two simple floating-point based routines (no explicit manipulation
of the internal double's format is performed).
*/
EightByteDouble nextRand_R0(void);
EightByteDouble nextRand_R1(void);
/*
Tests to see if a machine is big- or little-endian.
*/
int testEndian(void);
/*
A test routine to verify that mod 2^46 arithmetic is being performed
properly with this code on your machine.
*/
int testSeq46(nextRandType);
/*
A timing routine for the nextRand* routines.
*/
int testTime(nextRandType, skipNRandType);
/*
A routine to test how your compiler encodes floating-point
constants.
*/
int testCompiler46(nextRandType);
/*
A routine to test skipping N values in the series.
*/
int testNth(nextRandType, skipNRandType);
/*
A routine to test skipping N values in the series, for very large N.
*/
int testLargeNth(nextRandType, skipNRandType);
/*
This routine runs the whole battery of tests.
*/
int runTest(void);
/*
This routine demonstrates sample use of the generator.
*/
int SampleUse(void);
/*
This array will contain masks used in clearing the low bits of a
mantissa. (The values are computed by setMantissaMasks.)
*/
static int mantissaMasks[64] = {
0x00000000, 0xfff00000, 0x00000000, 0xfff80000,
0x00000000, 0xfffc0000, 0x00000000, 0xfffe0000,
0x00000000, 0xffff0000, 0x00000000, 0xffff8000,
0x00000000, 0xffffc000, 0x00000000, 0xffffe000,
0x00000000, 0xfffff000, 0x00000000, 0xfffff800,
0x00000000, 0xfffffc00, 0x00000000, 0xfffffe00,
0x00000000, 0xffffff00, 0x00000000, 0xffffff80,
0x00000000, 0xffffffc0, 0x00000000, 0xffffffe0,
0x00000000, 0xfffffff0, 0x00000000, 0xfffffff8,
0x00000000, 0xfffffffc, 0x00000000, 0xfffffffe,
0x00000000, 0xffffffff, 0x80000000, 0xffffffff,
0xc0000000, 0xffffffff, 0xe0000000, 0xffffffff,
0xf0000000, 0xffffffff, 0xf8000000, 0xffffffff,
0xfc000000, 0xffffffff, 0xfe000000, 0xffffffff,
0xff000000, 0xffffffff, 0xff800000, 0xffffffff,
0xffc00000, 0xffffffff, 0xffe00000, 0xffffffff};
#if 1
/*
This routine computes the masks used to clear out the desired bits
of a mantissa. This routine need not be called, as the values are
defined as part of the source.
*/
void setMantissaMasks(void);
void setMantissaMasks()
{
int i;
for (i=0; i<32; i++)
mantissaMasks[i*2+1] = 0xFFFFFFFF << ((i<20) ? (20 - i) : 0);
for (i=0; i<32; i++)
mantissaMasks[i*2] = (i>20) ? (0xFFFFFFFF << (52 - i)) : 0;
printf("static int mantissaMasks[64] = {\n");
for (i=0; i<63; i++)
{
printf("0x%08x, ", mantissaMasks[i]);
if (((i-3) % 4) == 0)
printf("\n");
}
printf("0x%08x};\n", mantissaMasks[63]);
}
#endif
static fp seed;
void setRandSeed(EightByteDouble s)
{
seed.u.d = s;
}
EightByteDouble getRandSeed()
{
return(seed.u.d);
}
#if PARTSROUTINE
/*
This routine divides the dividend by 2^power, storing the result
in the quotient and the residual in dividend_residual.
This is used only for debugging, as the macro (below) is used
"in-line" rather than calling a routine.
*/
void PartsR(fp *dividend_residual, fp *quotient, int power)
{
register int EXP_BASE = EXP(*dividend_residual) - power - BASE;
/*
If exp - power - base is less than 0, then the quotient is 0, and
the dividend_residual need not be changed (as the dividend is also
the residual).
*/
if (EXP_BASE < 0)
quotient->u.d = 0;
else
{
/*
Setup quotient as
floor(dividend_residual/(2^power)) * 2^power
*/
KEEP_N_PLUS_1_BITS(*dividend_residual, *quotient, EXP_BASE);
if (quotient->u.d != 0)
{
/*
Subtract the quotient*2^power from the dividend, resulting in the
residual.
*/
dividend_residual->u.d -= quotient->u.d;
/*
Compute the actual quotient by dividing by 2^power
*/
EXP(*quotient) = EXP_BASE + BASE;
}
}
}
#endif
/*
The macro PartsM, is identical to PartsR, above, but is used for
in-line expansion.
*/
#define PartsM(dividend_residual, quotient, POWER) \
{\
register int EXP_BASE = EXP(dividend_residual) - (POWER) - BASE;\
\
if (EXP_BASE < 0)\
quotient.u.d = 0;\
else\
{\
KEEP_N_PLUS_1_BITS(dividend_residual, quotient, EXP_BASE);\
if (quotient.u.d != 0)\
{\
dividend_residual.u.d -= quotient.u.d;\
EXP(quotient) = EXP_BASE + BASE;\
}\
}\
\
}
#if PARTSROUTINE
/*
This variant of PartsR (and the matching macro below) return the
integer quotient multiplied by 2^power, rather than the true
quotient. This is an optimization used by nextRandFinal. On some
platforms this yields a faster result.
*/
void PartsR2(fp *dividend_residual, fp *quotient, int power)
{
register int EXP_BASE = EXP(*dividend_residual) - power - BASE;
/*
If exp - power - base is less than 0, then the quotient is 0, and
the dividend_residual need not be changed (as the dividend is also
the residual).
*/
if (EXP_BASE < 0)
quotient->u.d = 0;
else
{
/*
Setup quotient as
floor(dividend_residual/(2^power)) * 2^power
*/
KEEP_N_PLUS_1_BITS(*dividend_residual, *quotient, EXP_BASE);
/*
Subtract the quotient*2^power from the dividend, resulting in the
residual.
*/
dividend_residual->u.d -= quotient->u.d;
}
}
#endif
/*
The macro PartsM2, is identical to PartsR2, above, but is used for
in-line expansion.
*/
#define PartsM2(dividend_residual, quotient, POWER) \
{\
register int EXP_BASE = EXP(dividend_residual) - (POWER) - BASE;\
\
if (EXP_BASE < 0)\
quotient.u.d = 0;\
else\
{\
KEEP_N_PLUS_1_BITS(dividend_residual, quotient, EXP_BASE);\
dividend_residual.u.d -= quotient.u.d;\
}\
\
}
/*
mod2N is a shortened version of the Parts code above, and only
computes A = A mod (2^power). This is used when the quotient is
not required.
*/
#define mod2N(A, POWER) \
{\
fp TEMP_MOD2N;\
register int EXP_BASE = EXP(A) - (POWER) - BASE;\
\
if (EXP_BASE >= 0)\
{\
KEEP_N_PLUS_1_BITS(A, TEMP_MOD2N, EXP_BASE);\
(A).u.d -= TEMP_MOD2N.u.d;\
}\
\
}
/*
mod2Nf is a variation on mod2N that computes A mod 2^POWER.
This variant also computes (A mod 2^POWER)/2^POWER, returning
a fraction, which is often computed following a call to mod2N.
HINT: this macro uses KEEP_N_PLUS_1_BITS_SMALL,
which assumes that the third argument is NEVER greater
than 20.
Also, when we divide AFrac (in the last line
of this macro) by modifying its exponent directly, we
know that AFrac will never be 0, so we need not check
its value first.
If this code is used in another setting where A may be a multiple
of 2^POWER, then we must add code to check if AFrac is 0 first,
and skip the division if it is equal to 0.
*/
#define mod2Nf(A, AFrac, POWER) \
{\
fp TEMP_MOD2N;\
register int exp_P = EXP(A) - ((POWER) + BASE);\
\
if (exp_P >= 0)\
{\
KEEP_N_PLUS_1_BITS_SMALL(A, TEMP_MOD2N, exp_P);\
(A).u.d -= TEMP_MOD2N.u.d;\
}\
(AFrac).u.d = (A).u.d;\
EXP(AFrac) = exp_P + BASE;\
\
}
/*
This is a generic routine for multiplying a and b and returning the
result mod 2^FIXED_POWER in c. It is assumed that a and b are both
less than 2^FIXED_POWER.
*/
void mulABmodP(fp a, fp b, fp *c)
{
fp aHi, bHi, t;
#if PARTSROUTINE
parts(&a, &aHi, FIXED_POWER_HALF);
parts(&b, &bHi, FIXED_POWER_HALF);
#else
parts(a, aHi, FIXED_POWER_HALF);
parts(b, bHi, FIXED_POWER_HALF);
#endif
t.u.d = aHi.u.d*b.u.d + a.u.d*bHi.u.d;
mod2N(t, FIXED_POWER_HALF);
EXP(t) += FIXED_POWER_HALF;
c->u.d = t.u.d + a.u.d*b.u.d;
mod2N(*c, FIXED_POWER);
}
/*
This routine squares a and returns the result
mod 2^FIXED_POWER in c.
*/
void squareModP(fp a, fp *c)
{
fp aHi, t;
#if PARTSROUTINE
parts(&a, &aHi, FIXED_POWER_HALF);
#else
parts(a, aHi, FIXED_POWER_HALF);
#endif
t.u.d = aHi.u.d*a.u.d;
mod2N(t, FIXED_POWER_HALF-1);
EXP(t) += FIXED_POWER_HALF+1;
c->u.d = t.u.d + a.u.d*a.u.d;
mod2N(*c, FIXED_POWER);
}
/*
Below are two routines that use only floating-point arithmetic
to perform the nextRand computation. The first is a straight port
from Fortran to C of the code distributed with the serial NAS
sample code. The second routine has a few optimizations.
*/
EightByteDouble nextRand_R0()
{
#define Power_half (1L< 0)
{
/*
Some compilers cache the value of temp, resulting in an old
version of temp being used.
*/
#if BUGGY_OPTIMIZER
temp.u.d = temp.u.d / ((1 << 16) *(double)(1 << 16));
#else
EXP(temp) -= 32;
#endif
n = temp.u.d;
Ptable = &table[32];
for (m=1; m <= n && m != 0; m <<= 1)
{
if (m & n)
mulABmodP(*Ptable, b, &b);
Ptable++;
}
}
return(b.u.d);
}
enum {
err_bad_base = 1,/* #define BASE is not set properly */
err_bad_size, /* struct size != double */
err_bad_struct, /* fp representation has unanticipated structure */
err_bad_times2, /* 1*2 does not increment exponent by 1 */
err_bad_seq, /* sequence is not being computed properly */
err_bad_comp, /* bad compiler bug (converting constant to fp) */
err_bad_nth, /* bad Nth seed bug */
err_bad_100knth, /* bad Nth seed bug (100k) */
err_bad_endian, /* Wrong endian setting (see #define BIGENDIAN) */
err_bad_large_n, /* the skip N seeds routine fails for very large N\n
(try setting BUGGY_OPTIMIZER to 1) */
err_bad_shift /* bad shift computation (compiler/machine bug) */
};
static char *errorStrings[] = {
"#define BASE is not set properly",
"struct size != double",
"fp representation has unanticipated structure",
"1*2 does not increment exponent by 1",
"sequence is not being computed properly",
"bad compiler bug (converting constant to fp)",
"bad Nth seed bug",
"bad Nth seed bug (100k)",
"Wrong endian setting (see #define BIGENDIAN)",
"the skip N seeds routine fails for very large N\n"
" (try setting BUGGY_OPTIMIZER to 1)",
"bad shift computation (compiler/machine bug)"};
/*
Test the first few numbers in a mod 2^46 series.
*/
int testSeq46(nextRandType nextRandf)
{
setRandSeed(INITIAL_SEED);
(*nextRandf)(); /* 55909509111989 */
if (LO_MANT(seed) != 0x82D40000 ||
HI_MANT(seed) != 0xCB65C9B8 ||
EXP(seed) != BASE+45) return (err_bad_seq);
(*nextRandf)(); /* 61155031930969 */
if (LO_MANT(seed) != 0x61640000 ||
HI_MANT(seed) != 0xDE7B0FD1 ||
EXP(seed) != BASE+45) return (err_bad_seq);
(*nextRandf)(); /* 45573031421645 */
if (LO_MANT(seed) != 0x9B340000 ||
HI_MANT(seed) != 0xA5CB3165 ||
EXP(seed) != BASE+45) return (err_bad_seq);
(*nextRandf)(); /* 55279057169489 */
if (LO_MANT(seed) != 0xB1440000 ||
HI_MANT(seed) != 0xC91AA243 ||
EXP(seed) != BASE+45) return (err_bad_seq);
(*nextRandf)(); /* 1250169187877 */
if (LO_MANT(seed) != 0x12800000 ||
HI_MANT(seed) != 0x9189F1F7 ||
EXP(seed) != BASE+40) return (err_bad_seq);
(*nextRandf)(); /* 27551538756233 */
if (LO_MANT(seed) != 0xD4480000 ||
HI_MANT(seed) != 0xC876BD71 ||
EXP(seed) != BASE+44) return (err_bad_seq);
(*nextRandf)(); /* 56099529399485 */
if (LO_MANT(seed) != 0x72F40000 ||
HI_MANT(seed) != 0xCC16C216 ||
EXP(seed) != BASE+45) return (err_bad_seq);
(*nextRandf)(); /* 28525870915841 */
if (LO_MANT(seed) != 0x08080000 ||
HI_MANT(seed) != 0xCF8D9339 ||
EXP(seed) != BASE+44) return (err_bad_seq);
(*nextRandf)(); /* 35461206354069 */
if (LO_MANT(seed) != 0xA2540000 ||
HI_MANT(seed) != 0x8101D26E ||
EXP(seed) != BASE+45) return (err_bad_seq);
return(0);
}
/*
Test the timing for calls to nextRandf.
*/
int testTime(nextRandType nextRandf, skipNRandType skipNRandf)
{
fp f;
FourByteInt i;
clock_t t;
t = clock();
f.u.d = INITIAL_SEED;
setRandSeed(f.u.d);
for (i=0; i<102400; i++)
(*nextRandf)();
if (LO_MANT(seed) != SEED_100K_LO ||
HI_MANT(seed) != SEED_100K_HI ||
EXP(seed) != BASE+SEED_100K_M) return (err_bad_100knth);
f.u.d = (*skipNRandf)(INITIAL_SEED, 102400);
if (LO_MANT(f) != SEED_100K_LO ||
HI_MANT(f) != SEED_100K_HI ||
EXP(f) != BASE+SEED_100K_M) return (err_bad_100knth);
t -= clock();
return(t);
}
/*
Test how a compiler converts fp constants and compare with
the mod 2^46 series.
*/
int testCompiler46(nextRandType nextRandf)
{
setRandSeed(INITIAL_SEED);
(*nextRandf)(); /* 55909509111989 */
if (seed.u.d != 5.5909509111989e+13) return (err_bad_comp);
(*nextRandf)(); /* 61155031930969 */
if (seed.u.d != 6.1155031930969e+13) return (err_bad_comp);
(*nextRandf)(); /* 45573031421645 */
if (seed.u.d != 4.5573031421645e+13) return (err_bad_comp);
(*nextRandf)(); /* 55279057169489 */
if (seed.u.d != 5.5279057169489e+13) return (err_bad_comp);
(*nextRandf)(); /* 1250169187877 */
if (seed.u.d != 1.250169187877e+12) return (err_bad_comp);
(*nextRandf)(); /* 27551538756233 */
if (seed.u.d != 2.7551538756233e+13) return (err_bad_comp);
(*nextRandf)(); /* 56099529399485 */
if (seed.u.d != 5.6099529399485e+13) return (err_bad_comp);
(*nextRandf)(); /* 28525870915841 */
if (seed.u.d != 2.8525870915841e+13) return (err_bad_comp);
(*nextRandf)(); /* 35461206354069 */
if (seed.u.d != 3.5461206354069e+13) return (err_bad_comp);
return(0);
}
int testNth(nextRandType nextRandf, skipNRandType skipNRandf)
{
EightByteDouble initialSeed, newSeed, skippedSeed;
int i, j, count;
/*
Test that skip(skip(2^30 - 1)) == skip(2^30)
*/
skippedSeed = (*skipNRandf)(INITIAL_SEED, 1L<<30);
newSeed = (*skipNRandf)((*skipNRandf)(INITIAL_SEED, (1L<<30)-1), 1);
if (skippedSeed != newSeed) return (err_bad_nth);
setRandSeed(INITIAL_SEED);
for (i=0; i< 23; i++)
{
count = 1001*(*nextRandf)();
initialSeed = getRandSeed();
for (j=1; j <= count; j++)
(*nextRandf)();
newSeed = getRandSeed();
skippedSeed = (*skipNRandf)(initialSeed, count);
if (skippedSeed != newSeed ) return(err_bad_nth);
}
return(0);
}
/*
This routine tests the skipNRand routine with values larger
than 2^32.
*/
int testLargeNth(nextRandType nextRandf, skipNRandType skipNRandf)
{
EightByteDouble initialSeed, newSeed, skippedSeed, largeN;
int i, j, count;
UFourByteInt largeInt = 0x80000000;
largeN = largeInt;
setRandSeed(INITIAL_SEED);
for (i=0; i<23; i++)
{
count = 101*(*nextRandf)()+101*(*nextRandf)();
initialSeed = getRandSeed();
for (j=0; j* 0) return(error);
t0 = -error;
error = testTime(nextRand_R1, skipNRandCalls0);
if (error > 0) return(error);
t1 = -error;
error = testTime(nextRandEarly, skipNRandCalls0);
if (error > 0) return(error);
t2 = -error;
error = testTime(nextRandSecond, skipNRandCalls0);
if (error > 0) return(error);
t3 = -error;
error = testTime(nextRandFinal, skipNRandCalls0);
if (error > 0) return(error);
t4 = -error;
printf("The time for testing nextRand_R0 was %d\n", t0);
printf("The time for testing nextRand_R1 was %d\n", t1);
printf("The time for testing nextRandEarly was %d\n", t2);
printf("The time for testing nextRandSecond was %d\n", t3);
printf("The time for testing nextRandFinal was %d\n", t4);
return (0);
}
int SampleUse()
{
int n;
double largeNum;
EightByteDouble sampleSeed;
/*
Set the initial seed and store it into a variable in order
to reset the generator's state at a later time.
*/
setRandSeed(INITIAL_SEED);
sampleSeed = getRandSeed();
/*
Note that nextRandFinal (and other members of its family) will
return a value between 0 and 1. (Neither 0 or 1 will ever
be returned.)
*/
/*
Get a random number between 0 and 999 (inclusive).
*/
n = 1000*nextRandFinal();
/*
Restore the seed.
*/
setRandSeed(sampleSeed);
/*
Get a random number between 0 and 1000 (inclusive).
*/
n = 1001*nextRandFinal();
/*
Get a random number in the range: (5, 6, 7, 8, 9, 10).
*/
n = 5 + 6*nextRandFinal();
/*
Skip ahead 100 numbers in the sequence.
*/
sampleSeed = skipNRandCalls0(getRandSeed(), 100);
setRandSeed(sampleSeed);
/*
Skip ahead 2^35+1 numbers in the sequence.
NOTE: Use caution when creating values larger than 2^32.
*/
largeNum = 1 + (1<<5) * ((double)(1L<<30));
sampleSeed = skipNRandCalls0(getRandSeed(), largeNum);
setRandSeed(sampleSeed);
n = 1001*nextRandFinal();
return(n);
}
#if STANDALONE
main()
{
int error;
/* setMantissaMasks(); */
error = runTest();
if (error)
{
printf(
"Bug: %s\n\n", errorStrings[error-1]);
printf(
"If this program reports a bug (with no suggested fix), the\n");
printf(
"first thing to do is to change the current selection for\n");
printf(
"KEEP_N_PLUS_1_BITS. There are multiple #define's for this\n");
printf(
"macro and some do not work with various compilers.\n\n");
printf(
"There are a few other #define's that may help out too, check\n");
printf(
"the code for helpful hints marked with \"HINT\".\n\n");
printf(
"If all else fails, contact the author: dantowitz@aol.com\n");
exit(1);
}
SampleUse();
return(0);
}
#endif
*