Page 214

C: Early Marsaglia Random Number Generator

   double rndm(long ss)
   //call with ss!=0 to seed, else with ss=0
     static int nBits = 32; // set to # bits used for integers
     static long m1 = 2^{nBits-2} + 2^{nBits-2}-1;
     static long m[17]; // 17 components must be initialized
     static long i, j;
     long seed, k, k0, k1, j0, j1, m2 = $2^{nBits/2};
     if( ss != 0 ) {
       if( ss < 0 ) ss = -ss;
       seed = ss % m1;
       if( seed % 2 == 0 ) seed = seed-1;
       k0 = 9069 % m2;  k1 = 9069/m2;
       j0 = seed % m2;  j1 = seed/m2;
       loop p = 0,1,...,16
         seed = j0*k0;
         j1 = ((seed/m2)+j0*k1+j1*k0)% (m2/2);
         j0 = seed % m2;
         m[p] = j0+m2*j1;
       end loop
       i = 4; j=16;
     }
     k=m[i]-m[j];
     if( k<0 ) k = k+m1;
     m[j] = k;
     i = i-1; if( i<0) i=16;
     j = j-1; if( j<0) i=16;
     return (double)k/m1;