TweetFollow Us on Twitter

URandomLib

Volume Number: 14 (1998)
Issue Number: 10
Column Tag: Tools Of The Trade

URandomLib: The Ultimate Macintosh Random-Number Generator

by Michael McLaughlin, McLean, VA

Include this class in your projects and never have to worry about random numbers again

The Value of Nothing

Try to think of nothing. It's difficult. Sensory data alone tend to bias our thoughts and the brain tries to perceive patterns in this stream even when there is nothing there.

Random numbers are the software analogue of nothing, the sound of no hands clapping. They are used primarily as input, either by themselves or in conjunction with other data.

The unique value of random input is that it is completely neutral. Patterns of any kind, discernable in the output, could not have come from such input and must, instead, be attributed to whatever additional systems are present. Typically, it is the behavior of these systems that is of interest and a random input stream is a way of exercising the software without telling it what to do.

Small wonder, then, that the generation of "random" numbers has always been, and continues to be, a perennial topic in computer science. Applications range from the trivial (e.g., games) to the deadly serious (e.g., Monte Carlo simulations of nuclear reactors). In the latter case, the quality of the random numbers is very important. This is one time when "rolling-your-own" is definitely not recommended.

Of course, any algorithm purporting to produce random numbers cannot really do so. At best, the output will be pseudo-random, meaning only that there are no detectable patterns in it. Tests for such patterns are an active area of research and can be quite sophisticated. Our goals here are more modest and we shall focus on creating random numbers, not testing them.

The utility class, URandomLib, that is described in this article is a complete pseudo-random number generator (PRNG), implemented as a library. URandomLib makes the creation of random numbers about as trivial as one could wish, while assuring unsurpassed quality and execution speed.

The speed comes from the fact the low-level function responsible for the random stream is coded in optimized assembly language. The quality of the output comes from having a world-class algorithm which produces numbers that are very random.

How Random Are They?

They are so random that you can use any of the individual bits just as you would the entire output value of ordinary generators. This is unusual and most PRNGs come with dire warnings against breaking up a random binary word into separate pieces. As we shall see, URandomLib does so with impunity and even uses this as an additional mechanism to decrease execution time.

All PRNGs generate new random numbers using the previous one(s) as input, but there are many different algorithms. The most common, by far, are the multiplicative congruential generators. With these algorithms, each random integer, X, comes from the formula

X[i+1] = (a*X[i]) % m

where a and m are (unsigned long) constants.

However, just any old a and m will not do. If you simply make them up, your random numbers will not be very random.

Randomness is one of the two necessary features of any PRNG. The other is a long period, the length of the random sequence before the numbers start repeating themselves. Speed is a third feature, not absolutely necessary but highly desirable.

When you pick inferior values for a and m, you can get bad results. Once upon a time, there was a famous PRNG known as RANDU. Almost everybody used it. RANDU was a multiplicative congruential generator with a = 65539 and m = 2147483648. The value of m (= 0x80000000) was chosen because it makes the modulo operation very easy, especially in assembly language where you can do whatever you like. The value of a (= 0x10003) was reportedly chosen because its binary representation has only three 1-bits, making multiplication unusually fast. Today, RANDU is used only as an example of how not to construct a PRNG. We shall see why later, when we compare it to URandomLib.

The generator algorithm in URandomLib is known as "Ultra." It is a strenuously tested compound generator. In this case, the output from the first generator is XORed with the output of an independent generator which, all by itself, is quite a good PRNG.

The first PRNG in Ultra is a subtract-with-borrow (SWB) generator which works as follows: [See Marsaglia and Zaman, in Further Reading, for details.]

Let b = 232 and m = b37 - b24 - 1, a prime number. If X[0] ... X[36] are 37 integers in the closed interval [0, b-1], not all zero or b-1, and c the carry (or borrow) bit from the previous operation, then the sequence constructed using the recursion

X[n] = (X[n-24] - X[n-37] - c) % b

has a period of m-1, about 10356. There is a lot more theory involved, as well as tricky implementation details, and it is not obvious that the sequence so generated will appear random, but it does. After passing through the second generator, the final output is even more random, and the period increased to about 10366.

URandomLib

The class URandomLib is not intended to be instantiated by the user. In fact, the library will not work properly if you declare objects of this class. Instead, by including URandomLib.cp in your project, and URandomLib.h in the modules that reference it, there will be a single, global object, PRNG. The constructor for PRNG will be called prior to main() and the destructor called after main() exits. Consequently, the library will behave like a system resource and its functionality will be available at all times.

There are 17 functions available in URandomLib (see Listings 1 and 2). The multiplicity of return types allows the generator to extract, from the random array, only the number of bytes actually necessary to produce the desired result. This minimizes the frequency of Refill() calls, which further increases the speed of URandomLib. The functionality, and random output, of this library may be summarized as follows:

Initialization
This call is necessary only if you wish to start with known seeds. The default constructor initializes PRNG automatically, with random seeds. Both seeds must be greater than zero, else random seeds will be used. Initialize() also calls SaveStart().
Save/Recall
In order to reproduce a sequence of random numbers exactly, it is necessary to restore the PRNG to a previous state. SaveStart() and RecallStart() perform this function. If a filename is passed with SaveStart, the state will also be saved to a file. The filename is an optional parameter to RecallStart().
Integer
There are seven integer formats available, ranging from UShort7() to ULong32().
Boolean
UBoolean() returns true or false, using up only one random bit in the process.
Uniform
There are four random uniform functions, two returning float precision and two double precision. (Usually, floats are cast to doubles.) Uniform_0_1() returns a U(0, 1) float; Uniform_m1_1() returns a U(-1, 1) float. In both cases, the return value has full precision no matter how small it is. Also, neither function ever returns zero or one. With the double-precision counterparts, DUniform_0_1() and DUniform_m1_1(), a zero value is an extremely remote possibility.
Normal
Normal(float mu, float sigma) returns a true normal (gaussian) variate, with mean = mu and standard deviation = sigma. Sigma must be greater than zero (not checked).
Expo
Expo(float lambda) returns an exponential variate, with mean = standard deviation = lambda. Lambda must be greater than zero (not checked).

Note that URandomLib usually returns floats, not doubles. This is done for speed (floats can fit into a register; doubles typically cannot). However, this is not much of a sacrifice since double-precision random quantities are rarely necessary. To get type double, the output of URandomLib can always be cast. For the same reason, the scale parameters of Normal and Expo are not checked.

Now it is time to see what we get for our money!

Pop Quiz

The program URandomLibTest (see Listing 3) exercises all of the functions of URandomLib, using known seeds. This provides a check for proper implementation. Most of this program was coded in C to illustrate that mixing C and C++ is straightforward.

In addition, a comparison with RANDU is carried out, testing the randomness of individual bits. This is done via CoinFlipTest, a simulation in which ten coins are flipped repeatedly in an attempt to reproduce the theoretical outcome, given by the tenth row of Pascal's Triangle, viz.,

1 10 45 120 210 252 210 120 45 10 1

The kth row of Pascal's Triangle gives the relative frequencies for the number of Heads [0-k] in a random trial using k coins. The sum along any row is 2k (here, 1024). Therefore, in this simulation, any integer multiple of 1024 trials will give integral expected frequencies, making this little quiz easy to grade.

The grade will be determined using the famous ChiSquare test. The ChiSquare statistic is computed as follows:

where o[k] and e[k] are the observed and expected frequencies for bin k, resp., and where the summation includes all frequency bins.

The nice thing about the ChiSquare statistic is that it is very easy to assess the difference between theory (expectation) and experiment. In this case, there are ten degrees-of-freedom, df, and the improbability of a given ChiSquare value is a known function of df. For instance, there is only a 5-percent chance of ChiSquare(10) > 18.3 if the results of this simulation are truly random. Additional critical points can be found in Listing 3.

Needless to say, URandomLib passes the CoinFlip test with flying colors whereas most other generators, including RANDU, do not. Check it out! It should be noted that this simulation is not a particularly difficult quiz for a PRNG. For examples of more stringent tests, read the classic discussion by Knuth (see Further Reading) and examine the Diehard test suite at http://stat.fsu.edu/~geo/diehard.html.

As indicated above, the development of PRNGs is a continuing area of research and URandomLib is clearly not the final word on the subject. Nevertheless, you will find it very hard to beat.

Listing 1: URandomLib.h

#pragma once

#ifndef __TYPES__   
#include <Types.h>   // to define Boolean
#endif

class URandomLib {
public:
   URandomLib();
    ~URandomLib() {};

   long ULong32();         // U[-2147483648, 2147483647]
   long ULong31();         // U[0, 2147483647]

   short UShort16();       // U[-32768, 32767]
   short UShort15();       // U[0,32767]

   short UShort8();        // (short) U[-128, 127]
   short UShort8u();       // (short) U[0, 255]
   short UShort7();        // (short) U[0, 127]

   Boolean UBoolean();     // true or false

   float Uniform_0_1();    // U(0,1) with >= 25-bit mantissa
   float Uniform_m1_1();   // U(-1,1), but excluding zero

   double DUniform_0_1();  // U[0,1) with <= 63-bit mantissa
   double DUniform_m1_1(); // U(-1,1) with <= 63-bit mantissa

   float Normal(float mu, float sigma);  // Normal(mean, std. dev. > 0)
   float Expo(float lambda);             // Exponential(lambda > 0)

   Boolean SaveStart(char *pathname = nil);
   Boolean RecallStart(char *pathname = nil);
 
   void      Initialize(unsigned long seed1 = 0,
                   unsigned long seed2 = 0);

private:
   void      Refill();   // low-level core routine

   struct {
      double               gauss;
      unsigned long   FSR[37], SWB[37], brw, seed1, seed2;
      long                  bits;
      short               byt, bit;
      char                  *ptr;
   }   Ultra_Remember;      // to restart PRNG from a known beginning

   double               Ultra_2n63, Ultra_2n31, Ultra_2n7,
                        Ultra_gauss;     // remaining gaussian variate
   unsigned long        Ultra_seed2;
   long                 Ultra_bits;      // bits for UBoolean
   short                Ultra_bit;       // # bits left in bits
};

static URandomLib   PRNG;

Listing 2: URandomLib.cp

#include <stdio.h>
#include <OSUtils.h>          // for GetDateTime()
#include <math.h>

#include "URandomLib.h"

unsigned long   Ultra_FSR[37],      // final random numbers
                     Ultra_SWB[37], // subtract-with-borrow output
                     Ultra_brw,     // either borrow(68K) or ~borrow(PPC)
                     Ultra_seed1;
short               Ultra_byt;      // # bytes left in FSR[37]
char                  *Ultra_ptr;   // running pointer to FSR[37]

Constructor, Destructor
URandomLib is initialized with random seeds, based on the system clock. There is a stub destructor.

URandomLib::URandomLib()
{
   Initialize();
}

URandomLib::~URandomLib() {};

Refill
This is the core of URandomLib. It refills Ultra_SWB[37] via a subtract-with-borrow PRNG, then superimposes a multiplicative congruential PRNG to produce Ultra_FSR[37], which supplies all of the random bytes.

#if defined(powerc)
asm void URandomLib::Refill()
{
      lwz      r3,Ultra_brw      // fetch global addresses from TOC
      lwz      r6,Ultra_SWB
      lwz      r4,0(r3)          // ~borrow
      la         r7,48(r6)       // &Ultra_SWB[12]
      sub      r5,r5,r5          // clear entire word
      mr         r8,r5           // counter
      li         r5,1
      sraw      r4,r4,r5         // restore XER|CA
      li         r8,24
      mtctr   r8
      la         r4,-4(r6)
UR1:   lwzu      r9,4(r7)
      lwz      r10,4(r4)
      subfe   r9,r10,r9          // r9 -= r10
      stwu      r9,4(r4)
      bdnz+   UR1
      mr         r7,r6           // &Ultra_SWB
      li         r8,13
      mtctr   r8
      la         r7,-4(r6)
UR2:   lwzu      r9,4(r7)
      lwz      r10,4(r4)
      subfe   r9,r10,r9          // r9 -= r10
      stwu      r9,4(r4)
      bdnz+   UR2
      lwz      r4,0(r3)          // ~borrow again
      addme   r4,r5              // r5 = 1
      neg      r4,r4
      stw      r4,0(r3)          // new ~borrow
      la         r6,-4(r6)       // &SWB[-1]
      lwz      r7,Ultra_FSR
      lwz      r5,Ultra_ptr
      lwz      r4,Ultra_seed1
      stw      r7,0(r5)          // reset running pointer to FSR
      la         r7,-4(r7)       // overlay congruential PRNG
      lis      r10,1             // r10 = 69069
      addi      r10,r10,3533
      lwz      r5,0(r4)          // Ultra_seed1
      li         r8,37
      mtctr   r8
UR3:   lwzu      r9,4(r6)        // SWB
      mullw   r5,r5,r10          // Ultra_seed1 *= 69069
      xor      r9,r9,r5
      stwu      r9,4(r7)            
      bdnz+   UR3
      stw      r5,0(r4)          // save Ultra_seed1 for next time
      lwz      r7,Ultra_byt
      li         r5,148          // 4*37 bytes
      sth      r5,0(r7)          // reinitialize
      blr
}
#else
asm void URandomLib::Refill()
{
      machine   68020

      MOVE.L   A2,-(SP)          // not scratch
      LEA      Ultra_SWB,A2      // &Ultra_SWB[0]
      LEA      52(A2),A1         // &Ultra_SWB[13]
      MOVEQ   #0,D0              // restore extend bit
      SUB.L   Ultra_brw,D0
      MOVEQ   #23,D2             // 24 of these
UR1:   MOVE.L   (A1)+,D0
      MOVE.L   (A2),D1
      SUBX.L   D1,D0
      MOVE.L   D0,(A2)+
      DBRA      D2,UR1
      LEA      Ultra_SWB,A1
      MOVEQ   #12,D2             // 13 of these
UR2:   MOVE.L   (A1)+,D0
      MOVE.L   (A2),D1
      SUBX.L   D1,D0             // subtract-with-borrow
      MOVE.L   D0,(A2)+
      DBRA      D2,UR2
      MOVEQ   #0,D0
      MOVE.L   D0,D1
      ADDX      D1,D0            // get borrow bit
      MOVE.L   D0,Ultra_brw      //   and save it
      LEA      Ultra_SWB,A1
      LEA      Ultra_FSR,A2
      MOVE.L   A2,Ultra_ptr      // reinitialize running pointer
      MOVE.L   Ultra_seed1,D0
      MOVE.L   #69069,D1         // overlay congruential PRNG
      MOVEQ   #36,D2             // 37 of these
UR3:   MOVE.L   (A1)+,(A2)
      MULU.L   D1,D0
      EOR.L   D0,(A2)+
      DBRA      D2,UR3
      MOVE.L   D0,Ultra_seed1    // save global for next time
      MOVE      #148,Ultra_byt   // 4*37 bytes left
      MOVE.L   (SP)+,A2          // restore A2
      RTS
}
#endif

ULong32
ULong32() returns a four-byte integer, ~U[-2147483648, 2147483647]. It may, of course, be cast to unsigned long.

long URandomLib::ULong32()
{
   register long   result;
   
   if (Ultra_byt < 4) Refill();
   result = *((long *) Ultra_ptr);
   Ultra_ptr += 4; Ultra_byt -= 4;
   return result;
}

ULong31
ULong31() returns a four-byte integer, ~U[0, 2147483647].

long URandomLib::ULong31()
{
   register long   result;
   
   if (Ultra_byt < 4) Refill();
   result = *((long *) Ultra_ptr);
   Ultra_ptr += 4; Ultra_byt -= 4;
   return result & 0x7FFFFFFF;
}

UShort16
UShort16() returns a two-byte integer, ~U[-32768, 32767].

short URandomLib::UShort16()
{
   register short   result;
   
   if (Ultra_byt < 2) Refill();
   result = *((short *) Ultra_ptr);
   Ultra_ptr += 2; Ultra_byt -= 2;
   return result;
}

UShort15
UShort15() returns a two-byte integer, ~U[0, 32767].

short URandomLib::UShort15()
{
   register short   result;
   
   if (Ultra_byt < 2) Refill();
   result = *((short *) Ultra_ptr);
   Ultra_ptr += 2; Ultra_byt -= 2;
   return result & 0x7FFF;
}

UShort8
UShort8() returns a two-byte integer, ~U[-128, 127]. It gets a random byte and casts it to short. This operation extends the sign bit. Consequently, you may NOT cast this function to unsigned short/long (see UShort8u() below).

short URandomLib::UShort8()
{
   register short   result;
   
   if (Ultra_byt < 1) Refill();
   result = (short) *Ultra_ptr;
   Ultra_ptr += 1; Ultra_byt -= 1;
   return result;
}

UShort8u
UShort8u() returns a two-byte integer, ~U[0, 255]. It proceeds as in UShort8() but clears the high byte instead of extending the sign bit.

short URandomLib::UShort8u()
{
   register short   result;
   
   if (Ultra_byt < 1) Refill();
   result = (short) *Ultra_ptr;
   Ultra_ptr += 1; Ultra_byt -= 1;
   return result & 0xFF;
}

UShort7
UShort7() returns a two-byte integer, ~U[0, 127].

short URandomLib::UShort7()
{
   register short   result;
   
   if (Ultra_byt < 1) Refill();
   result = (short) (*Ultra_ptr & 0x7F);
   Ultra_ptr += 1; Ultra_byt -= 1;
   return result;
}

UBoolean
UBoolean() returns true or false. It calls ULong32() and returns the bits one at a time.

Boolean URandomLib::UBoolean()
{
   register Boolean   result;
   
   if (Ultra_bit <= 0) {
      Ultra_bits = ULong32();
      Ultra_bit = 32;
   }
   result = (Ultra_bits < 0) ? true : false;
   Ultra_bits += Ultra_bits;   // shift left by one
   &#151;Ultra_bit;
   return result;
}

Uniform_0_1
Uniform_0_1() returns a four-byte float, ~U(0, 1), with >= 25 bits of precision. This precision is achieved by continually testing the leading byte, b, of the mantissa. If b == 0, it is replaced with a new random byte and the decimal point readjusted. This simultaneously ensures that Uniform_0_1() never returns zero.

float URandomLib::Uniform_0_1()
{
   register double      fac = Ultra_2n31;
   register long      along;
   register short      extra;
   
   along = ULong31();
   if (along >= 0x01000000) return (float)(fac*along);
   for (extra=0;!extra;) {      // will not be an infinite loop
       extra = UShort7();
       fac *= Ultra_2n7;
   }
   along |= (((long)extra) << 24);
   return (float)(fac*along);
}

Uniform_m1_1
Uniform_m1_1() returns a four-byte float, ~U(-1, 1), with the same features as described above for Uniform_0_1().

float URandomLib::Uniform_m1_1()
{
   register double    fac = Ultra_2n31;
   register long      along, limit = 0x01000000;
   register short     extra;
   
   if ((along = ULong32()) >= limit)
      return (float)(fac*along);
   else if (-along >= limit)
      return (float)(fac*along);
   for (extra=0;!extra;) {
       extra = UShort7();
       fac *= Ultra_2n7;
   }
   if (along >= 0) {
      along |= (((long)extra) << 24);
      return (float)(fac*along);
   }
   along = -along;
   along |= (((long)extra) << 24);
   return (float)(-fac*along);
}

DUniform_0_1, DUniform_m1_1
DUniform_0_1() and DUniform_m1_1() return double-precision U[0,1) and U(-1,1). In both cases, zero IS a remote possibility. These functions are intended for those occasions when seven significant figures are not enough. If you need TYPE double, but not double PRECISION, then it is much faster to use Uniform_0_1() or Uniform_m1_1() and cast - implicitly or explicitly.

double URandomLib::DUniform_0_1()
{
   return ULong31()*Ultra_2n31 +
            ((unsigned long) ULong32())*Ultra_2n63;
}

double URandomLib::DUniform_m1_1()
{
   return ULong32()*Ultra_2n31 +
            ((unsigned long) ULong32())*Ultra_2n63;
}

Normal
Normal() returns a four-byte float, ~Normal(mu, sigma), where mu and sigma are the mean and standard deviation, resp., of the parent population. The normal variates returned are exact, not approximate. Normal() uses Uniform_m1_1() so there is no possibility of a result exactly equal to mu. Note that mu and sigma must also be floats, not doubles.

float URandomLib::Normal(float mu, float sigma)
{
   register double      fac, rsq, v1, v2;

   if ((v1 = Ultra_gauss) != 0.0) {      // Is there one left?
      Ultra_gauss = 0.0;
      return (float)(sigma*v1 + mu);
   }
   do {
      v1 = Uniform_m1_1();
      v2 = Uniform_m1_1();
      rsq = v1*v1 + v2*v2;
   } while (rsq >= 1.0);
   fac = sqrt(-2.0*log(rsq)/rsq);
   Ultra_gauss = fac*v2;                 // Save the first N(0,1) as double
   return (float)(sigma*fac*v1 + mu);    // and return the second
}

Expo
Expo() returns a four-byte float, ~Exponential(lambda). The parameter, lambda, is both the mean and standard deviation of the parent population. It must be a float greater than zero.

float URandomLib::Expo(float lambda)
{
   return (float)(-lambda*log(Uniform_0_1()));
}

SaveStart, RecallStart
SaveStart() and RecallStart() save and restore, resp., the complete state of URandomLib. Call SaveStart() at the point where it may be necessary to recall a sequence of random numbers exactly. To recover the sequence later, call RecallStart(). To terminate a program and still recover a random sequence, save Ultra_Remember to a file and read it back upon restart.

Boolean URandomLib::SaveStart(char *pathname)
{
   Ultra_Remember.gauss = Ultra_gauss;
   Ultra_Remember.bits = Ultra_bits;
   Ultra_Remember.seed1 = Ultra_seed1;
   Ultra_Remember.seed2 = Ultra_seed2;
   Ultra_Remember.brw = Ultra_brw;
   Ultra_Remember.byt = Ultra_byt;
   Ultra_Remember.bit = Ultra_bit;
   Ultra_Remember.ptr = Ultra_ptr;
   for (int i = 0;i < 37;i++) {
      Ultra_Remember.FSR[i] = Ultra_FSR[i];
      Ultra_Remember.SWB[i] = Ultra_SWB[i];
   }
   
   if (pathname != nil) {
      FILE   *outfile;
      if ((outfile = fopen(pathname, "w")) != nil) {
         fwrite((void *) &Ultra_Remember,
                  sizeof(Ultra_Remember), 1L, outfile);
         fclose(outfile);
      }
      else return false;
   }
   return true;
}

Boolean URandomLib::RecallStart(char *pathname)
{
   if (pathname != nil) {
      FILE   *infile;
      if ((infile = fopen(pathname, "r")) != nil) {
         fread((void *) &Ultra_Remember,
                  sizeof(Ultra_Remember), 1L, infile);
         fclose(infile);
      }
      else return false;
   }

   Ultra_gauss = Ultra_Remember.gauss;
   Ultra_bits = Ultra_Remember.bits;
   Ultra_seed1 = Ultra_Remember.seed1;
   Ultra_seed2 = Ultra_Remember.seed2;
   Ultra_brw = Ultra_Remember.brw;
   Ultra_byt = Ultra_Remember.byt;
   Ultra_bit = Ultra_Remember.bit;
   Ultra_ptr = Ultra_Remember.ptr;
   for (int i = 0;i < 37;i++) {
      Ultra_FSR[i] = Ultra_Remember.FSR[i];
      Ultra_SWB[i] = Ultra_Remember.SWB[i];
   }

   return true;
}

Initialize
Initialize() computes a few global constants, initializes others, and fills in the initial Ultra_SWB array using the supplied seeds. It terminates by calling SaveStart() so that you may recover the whole sequence of random numbers by calling RecallStart().

void URandomLib::Initialize(unsigned long seed1,
                            unsigned long seed2)
{
#if defined(powerc)
#define   ULTRABRW      0xFFFFFFFF
#else
#define   ULTRABRW      0x00000000
#endif

   unsigned long   tempbits, ul, upper, lower;
   
   if ((seed1 == 0) || (seed2 == 0)) {   // random initialization
      ::GetDateTime(&seed1);
      upper = (seed1 & 0xFFFF0000) >> 16;
      lower = seed1 & 0xFFFF;
      seed2 = upper*lower;               // might overflow
   }
   Ultra_seed1 = seed1; Ultra_seed2 = seed2;
    for (int i = 0;i < 37;i++) {
      tempbits = 0;
      for (int j = 32;j > 0;j&#151;) {
         Ultra_seed1 *= 69069;
         Ultra_seed2 ^= (Ultra_seed2 >> 15);
         Ultra_seed2 ^= (Ultra_seed2 << 17);
         ul = Ultra_seed1 ^ Ultra_seed2;
         tempbits = (tempbits >> 1) | (0x80000000 & ul);
      }
      Ultra_SWB[i] = tempbits;
   }
   Ultra_2n31 = ((2.0/65536)/65536);
   Ultra_2n63 = 0.5*Ultra_2n31*Ultra_2n31;
   Ultra_2n7 = 1.0/128;
   Ultra_gauss = 0.0;
   Ultra_byt = Ultra_bit = 0;
   Ultra_brw = ULTRABRW;                 // no borrow yet
   SaveStart();
}

Listing 3: URandomLibTest.cp

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "URandomLib.h"

/* Prototypes */
Boolean RANDU_Boolean();
void CoinFlipTest(int rpt, Boolean URLib);
double ChiSquare(long result[], int df);
double ExerciseAll();
void main();
/**************/

long      RANDU_Seed, Expectation[11],
         Theory[11] = {1,10,45,120,210,252,210,120,45,10,1};

CoinFlipTest
CoinFlipTest () attempts to reproduce an integer multiple (rpt) of the tenth row of Pascal's Triangle by flipping ten coins at a time.

void CoinFlipTest(int rpt, Boolean URLib)
{
   double   ans;
   long      i, PascalRow10[11];
   int       coin, Heads;
   
   static double crit[10] = 
   {3.94,16.0,18.3,23.2,29.6,35.6,41.3,46.9,52.3,57.7};
   static double conf[10] = 
   {0,90,95,99,99.9,99.99,99.999,99.9999,99.99999,99.999999};

   for (i = 0;i <= 10;i++)
      PascalRow10[i] = 0;
   if (URLib) {      // use URandomLib
      for (i = 1;i <= rpt*1024;i++) {
         Heads = 0;
         for (coin = 1;coin <= 10;coin++)
            if (PRNG.UBoolean()) ++Heads;
         ++PascalRow10[Heads];
      }
   }
   else {               // use RANDU
      for (i = 1;i <= rpt*1024;i++) {
         Heads = 0;
         for (coin = 1;coin <= 10;coin++)
            if (RANDU_Boolean()) ++Heads;
         ++PascalRow10[Heads];
      }
   }
   for (i = 0;i <= 10;i++)
      printf("%ld ", PascalRow10[i]);
   printf("\n\n");
   ans = ChiSquare(PascalRow10, 10);
   printf("ChiSquare = %f ==> ", ans);
   if (ans < crit[0])
      printf("Result is suspiciously good!\n\n");
   else if (ans > crit[1]) {
      int k;
      for (k = 1;(k <= 8) && (ans > crit[k+1]);) ++k;
      printf("Randomness is rejected with more than %f%% 
                  confidence.\n\n", conf[k]);
   }
   else printf("Randomness is accepted.\n\n");
}

ChiSquare
Compute the ChiSquare statistic for df degrees-of-freedom. The expected value = df.

double ChiSquare(long result[], int df)
{
   double   diff, chisq = 0.0;

   for (int i = 0;i <= df;i++) {
      diff = result[i] - Expectation[i];
      chisq += (diff*diff)/Expectation[i];
   }
   return chisq;
}

RANDU_Boolean
RANDU_Boolean() gets bits in much the same fashion as URandomLib.

Boolean RANDU_Boolean()
{
   Boolean   result;
   static unsigned long   a = 65539,      // RANDU constants
                                 m = 2147483648;
   static long theBits;
   static int bits_left = 0;
   
   if (bits_left <= 0) {
      theBits = RANDU_Seed =
                     (a*RANDU_Seed) % m;  // RANDU
      theBits += theBits;                 // initial sign bit always zero
      bits_left = 31;
   }
   result = (theBits < 0) ? true : false;
   theBits += theBits;                    // shift left by one
   &#151;bits_left;
   return result;
}

ExerciseAll
ExerciseAll () tests all of the functions in URandomLib.

double ExerciseAll()
{
   double   total = 0.0;
   float      mean, sigma;
   short      k;
   
   for (long i = 0;i < 50000;i++) {
      k = PRNG.UShort7() & 15;
      switch (k) {
         case 0:
            total += (double)PRNG.ULong32();
            break;
         case 1:
            total += (double)PRNG.ULong31();
            break;
         case 2:
            total -= (double)PRNG.ULong31();
            break;
         case 3:
            total += (double)PRNG.UShort16();
            break;
         case 4:
            total += (double)PRNG.UShort15();
            break;
         case 5:
            total -= (double)PRNG.UShort15();
            break;
         case 6:
            total += (double)PRNG.UShort8();
            break;
         case 7:
            total += (double)PRNG.UShort8u();
            break;
         case 8:
            total += (double)PRNG.UShort7();
            break;
         case 9:
            total += (double)PRNG.UBoolean();
            break;
         case 10:
            total += (double)PRNG.Uniform_0_1();
            break;
         case 11:
            total += (double)PRNG.Uniform_m1_1();
            break;
         case 12:
            total += (double)PRNG.DUniform_0_1();
            break;
         case 13:
            total += (double)PRNG.DUniform_m1_1();
            break;
         case 14:
            mean = PRNG.Uniform_m1_1();
            sigma = PRNG.Uniform_0_1();
            total += (double)PRNG.Normal(mean, sigma);
            break;
         case 15:
            total += (double)PRNG.Expo(PRNG.Uniform_0_1());
      }
   }
   return total;
}

main
Carry out CoinFlipTest and ExerciseAll.

void main()
{
   int   Nrepeats;
   
   // initialize RANDU
   RANDU_Seed = PRNG.ULong32();      // PRNG is automatically initialized
   // test individual "random" bits
   printf("Coin-flip test:\n\n");
   printf("Enter the number of repetitions.\n");
   scanf("%d", &Nrepeats);
   printf("\n");
   printf("Expected frequencies:\n");
   for (int i = 0;i <= 10;i++) {
      Expectation[i] = Nrepeats*Theory[i];
      printf("%ld ", Expectation[i]);
   }
   printf("\n\n");
   printf("Using URandomLib...\n");
   CoinFlipTest(Nrepeats, true);      // use URandomLib
   printf("Using RANDU...\n");
   CoinFlipTest(Nrepeats, false);     // use RANDU
   // test all of the functions in URandomLib
   printf("Exercise all functions: 
               (you should get 1.381345e+11, twice)\n\n");
   PRNG.Initialize(12345678, 87654321);
   PRNG.SaveStart("UltraTemp.dat");   // save initial state to file
   printf("%e\n", ExerciseAll());
   PRNG.RecallStart("UltraTemp.dat"); // initial state from file
   printf("%e\n", ExerciseAll());
}

Bibliography and References

  • Marsaglia, George and Arif Zaman. "A New Class of Random Number Generators", Annals of Applied Probability, vol. 1 No. 3 (1991), pp. 462-480.
  • Knuth, Donald E. The Art of Computer Programming, 2nd ed., vol. 2, Chap. 3, Addison-Wesley, 1981.

Michael McLaughlin, mpmcl@mitre.org, a former chemistry professor and Peace Corps volunteer, currently does R&D for future Air Traffic Control systems. He has been programming computers since 1965 but has long since forsaken Fortran, PLI, and Lisp in favor of C++ and assembly.

 

Community Search:
MacTech Search:

Software Updates via MacUpdate

Microsoft Remote Desktop 8.0.19 - Connec...
With Microsoft Remote Desktop, you can connect to a remote PC and your work resources from almost anywhere. Experience the power of Windows with RemoteFX in a Remote Desktop client designed to help... Read more
OmniGraffle 6.3 - Create diagrams, flow...
OmniGraffle helps you draw beautiful diagrams, family trees, flow charts, org charts, layouts, and (mathematically speaking) any other directed or non-directed graphs. We've had people use Graffle to... Read more
PDFKey Pro 4.3.2 - Edit and print passwo...
PDFKey Pro can unlock PDF documents protected for printing and copying when you've forgotten your password. It can now also protect your PDF files with a password to prevent unauthorized access and/... Read more
Ableton Live 9.2.2 - Record music using...
Ableton Live lets you create and record music on your Mac. Use digital instruments, pre-recorded sounds, and sampled loops to arrange, produce, and perform your music like never before. Ableton Live... Read more
Macs Fan Control 1.3.1.0 - Monitor and c...
Macs Fan Control allows you to monitor and control almost any aspect of your computer's fans, with support for controlling fan speed, temperature sensors pane, menu-bar icon, and autostart with... Read more
NetShade 6.3.1 - Browse privately using...
NetShade is an anonymous proxy and VPN app+service for Mac. Unblock your Internet through NetShade's high-speed proxy and VPN servers spanning seven countries. NetShade masks your IP address as you... Read more
Dragon Dictate 4.0.7 - Premium voice-rec...
With Dragon Dictate speech recognition software, you can use your voice to create and edit text or interact with your favorite Mac applications. Far more than just speech-to-text, Dragon Dictate lets... Read more
Persecond 1.0.2 - Timelapse video made e...
Persecond is the easy, fun way to create a beautiful timelapse video. Import an image sequence from any camera, trim the length of your video, adjust the speed and playback direction, and you’re done... Read more
GIMP 2.8.14p2 - Powerful, free image edi...
GIMP is a multi-platform photo manipulation tool. GIMP is an acronym for GNU Image Manipulation Program. The GIMP is suitable for a variety of image manipulation tasks, including photo retouching,... Read more
Sandvox 2.10.2 - Easily build eye-catchi...
Sandvox is for Mac users who want to create a professional looking website quickly and easily. With Sandvox, you don't need to be a Web genius to build a stylish, feature-rich, standards-compliant... Read more

ReBoard: Revolutionary Keyboard (Utilit...
ReBoard: Revolutionary Keyboard 1.0 Device: iOS Universal Category: Utilities Price: $1.99, Version: 1.0 (iTunes) Description: Do everything within the keyboard without switching apps! If you are in WhatsApp, how do you schedule a... | Read more »
Tiny Empire (Games)
Tiny Empire 1.1.3 Device: iOS Universal Category: Games Price: $2.99, Version: 1.1.3 (iTunes) Description: Launch cannonballs and blow tiny orcs into thousands of pieces in this intuitive fantasy-themed puzzle shooter! Embark on an... | Read more »
Astropad Mini (Productivity)
Astropad Mini 1.0 Device: iOS iPhone Category: Productivity Price: $4.99, Version: 1.0 (iTunes) Description: *** 50% off introductory price! ​*** Get the high-end experience of a Wacom tablet at a fraction of the price with Astropad... | Read more »
Emo Chorus (Music)
Emo Chorus 1.0.0 Device: iOS Universal Category: Music Price: $1.99, Version: 1.0.0 (iTunes) Description: Realistic Choir simulator ranging from simple Chorus emulation to full ensemble Choir with 128 members. ### introductory offer... | Read more »
Forest Spirit (Games)
Forest Spirit 1.0.5 Device: iOS Universal Category: Games Price: $2.99, Version: 1.0.5 (iTunes) Description: | Read more »
Ski Safari 2 (Games)
Ski Safari 2 1.0 Device: iOS Universal Category: Games Price: $1.99, Version: 1.0 (iTunes) Description: The world's most fantastical, fun, family-friendly skiing game is back and better than ever! Play as Sven's sister Evana, share... | Read more »
Lara Croft GO (Games)
Lara Croft GO 1.0.47768 Device: iOS Universal Category: Games Price: $4.99, Version: 1.0.47768 (iTunes) Description: Lara Croft GO is a turn based puzzle-adventure set in a long-forgotten world. Explore the ruins of an ancient... | Read more »
Whispering Willows (Games)
Whispering Willows 1.23 Device: iOS Universal Category: Games Price: $4.99, Version: 1.23 (iTunes) Description: **LAUNCH SALE 50% OFF** - Whispering Willows is on sale for 50% off ($4.99) until September 9th. | Read more »
Calvino Noir (Games)
Calvino Noir 1.1 Device: iOS iPhone Category: Games Price: $3.99, Version: 1.1 (iTunes) Description: The film noir stealth game. Calvino Noir is the exploratory, sneaking adventure through the 1930s European criminal underworld.... | Read more »
Angel Sword (Games)
Angel Sword 1.0 Device: iOS Universal Category: Games Price: $6.99, Version: 1.0 (iTunes) Description: Prepare to adventure in the most epic full scale multiplayer 3D RPG for mobile! Experience amazing detailed graphics in full HD.... | Read more »

Price Scanner via MacPrices.net

Apple and Cisco Partner to Deliver Fast-Lane...
Apple and Cisco have announced a partnership to create a “fast lane” for iOS business users by optimizing Cisco networks for iOS devices and apps. The alliance integrates iPhone with Cisco enterprise... Read more
Apple offering refurbished 2015 13-inch Retin...
The Apple Store is offering Apple Certified Refurbished 2015 13″ Retina MacBook Pros for up to $270 (15%) off the cost of new models. An Apple one-year warranty is included with each model, and... Read more
Apple refurbished 2015 MacBook Airs available...
The Apple Store has Apple Certified Refurbished 2015 11″ and 13″ MacBook Airs (the latest models), available for up to $180 off the cost of new models. An Apple one-year warranty is included with... Read more
21-inch iMacs on sale for up to $120 off MSRP
B&H Photo has 21″ iMacs on sale for up to $120 off MSRP including free shipping plus NY sales tax only: - 21″ 1.4GHz iMac: $999.99 $100 off - 21″ 2.7GHz iMac: $1199.99 $100 off - 21″ 2.9GHz iMac... Read more
5K iMacs on sale for up to $150 off MSRP, fre...
B&H Photo has the 27″ 3.3GHz 5K iMac on sale for $1899.99 including free shipping plus NY tax only. Their price is $100 off MSRP. They have the 27″ 3.5GHz 5K iMac on sale for $2149.99 $2199.99, $... Read more
1.4GHz Mac mini, refurbished, available for $...
The Apple Store has Apple Certified Refurbished 1.4GHz Mac minis available for $419. Apple’s one-year warranty is included, and shipping is free. Their price is $80 off MSRP, and it’s the lowest... Read more
iPad Air 2 on sale for up to $100 off MSRP
Best Buy has iPad Air 2s on sale for up to $100 off MSRP on their online store for a limited time. Choose free shipping or free local store pickup (if available). Sale prices available for online... Read more
MacBook Airs on sale for $100 off MSRP
Best Buy has MacBook Airs on sale for $100 off MSRP on their online store. Choose free shipping or free local store pickup (if available). Sale prices for online orders only, in-store prices may vary... Read more
Big Grips Lift Handle For iPad Air and iPad A...
KEM Ventures, Inc. which pioneered the extra-large, super-protective iPad case market with the introduction of Big Grips Frame and Stand in 2011, is launching Big Grips Lift featuring a new super-... Read more
Samsung Launches Galaxy Tab S2, Its Most Powe...
Samsung Electronics America, Inc. has announced the U.S. release of the Galaxy Tab S2, its thinnest, lightest, ultra-fast tablet. Blending form and function, elegant design and multitasking power,... Read more

Jobs Board

*Apple* Evangelist - JAMF Software (United S...
The Apple Evangelist is responsible for building and cultivating strategic relationships with Apple 's small and mid-market business development field teams. This Read more
*Apple* Retail - Multiple Positions (US) - A...
Sales Specialist - Retail Customer Service and Sales Transform Apple Store visitors into loyal Apple customers. When customers enter the store, you're also the Read more
*Apple* Desktop Analyst - KDS Staffing (Unit...
…field and consistent professional recruiting achievement. Job Description: Title: Apple Desktop AnalystPosition Type: Full-time PermanentLocation: White Plains, NYHot Read more
*Apple* Systems Engineer (Mclean, VA and NYC...
Title: Apple Systems Engineer (Mclean, VA and NYC) Location: United States-New York-New York-200 Park Ave (22005) Other Locations: United States-Virginia-Vienna-Towers Read more
*Apple* Systems Engineer (Mclean, VA and NYC...
…Assist in providing strategic direction and technical leadership within the Apple portfolio, including desktops, laptops, and printing environment. This person will Read more
All contents are Copyright 1984-2011 by Xplain Corporation. All rights reserved. Theme designed by Icreon.