home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / unix / volume28 / mrandom-3.0 / part01 / src / ultra.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-05-06  |  5.0 KB  |  172 lines
  1. /* ultra.c 3.1 5/28/93 */
  2. /*
  3. FSU - ULTRA    The greatest random number generator that ever was
  4.         or ever will be.  Way beyond Super-Duper.
  5.         (Just kidding, but we think its a good one.)
  6.  
  7. Authors:    Arif Zaman (arif@stat.fsu.edu) and
  8.         George Marsaglia (geo@stat.fsu.edu).
  9.  
  10. Date:        27 May 1992
  11.  
  12. Version:    1.05
  13.  
  14. Copyright:    To obtain permission to incorporate this program into
  15.         any commercial product, please contact the authors at
  16.         the e-mail address given above or at
  17.  
  18.         Department of Statistics and
  19.         Supercomputer Computations Research Institute
  20.         Florida State University
  21.         Tallahassee, FL 32306.
  22.  
  23. See Also:    README        for a brief description
  24.         ULTRA.DOC    for a detailed description
  25.  
  26. -----------------------------------------------------------------------
  27.                 Modified by Robert Plotkin for use with mrandom, 4/93
  28. */
  29. /*
  30.    File: ULTRA.C
  31.  
  32.    This is the ULTRA random number generator written entirely in C.
  33.  
  34.    This may serve as a model for an assembler version of this routine.
  35.    The programmer should avoid simply duplicating and instead use the
  36.    usual assembler features to increase the speed of this routine.
  37.  
  38.    Especially the subroutine SWB should be replaced by the one
  39.    machine instruction (usually called subtract-with-borrow) that
  40.    is available in almost every hardware.
  41.  
  42.    For people not familiar with 8086 assembler, it may help to
  43.    consult this when reading the assembler code. This program should
  44.    be a dropin replacement for the assembler versions, but is about
  45.    half as fast.
  46. */
  47.  
  48. #include "ultra.h"
  49.  
  50. #define N  37           /* size of table        */
  51. #define N2 24           /* The shorter lag      */
  52.  
  53. #define swb32   RNGstate
  54. #define swbseed (RNGstate+N)
  55. #define swb32p  ((long *)RNGstate[N+N])
  56. #define swb32n  RNGstate[N+N+1]
  57. #define flags   ((char)RNGstate[N+N+2]) /* The carry flag for
  58.                        the SWB generator    */
  59. #define congx   ((long unsigned)RNGstate[N+N+3]) /* Seed */
  60.  
  61. /* SWB is the subtract-with-borrow operation which should be one line
  62.    in assembler code. This should be done by using the hardware s-w-b
  63.    operation in the SWBfill routine.
  64.  
  65.    What has been done here is to look at the msb of x, y and z=x-y-c.
  66.    Using these three bits, one can determine if a borrow bit is needed
  67.    or not according to the following table:
  68.  
  69.     msbz=0  msby=0  msby=1          msbz=1  msby=0  msby=1
  70.  
  71.     msbx=0  0       1               msbx=0  1       1
  72.     msbx=1  0       0               msbx=1  0       1
  73.  
  74.    PS: note that the definition is very carefully written because the
  75.    calls to SWB have y and z as the same memory location, so y must
  76.    be tested before z is assigned a value.
  77. */
  78. #define SWB(c,x,y,z) \
  79. c = (y<0) ? (((z=x-y-c) < 0) || (x>=0)) : (((z=x-y-c) < 0) && (x>=0))
  80.  
  81. /*
  82.   The first two lines of this macro are the heart of the system and
  83.   should be written is assembler to be as fast as possible. It may even
  84.   make sense to unravel the loop and simply write 37 consecutive SWB
  85. operations!  */
  86. void SWBfill(x,rng)
  87. long *x;
  88. RNGdata *rng;
  89. { short i;
  90.  
  91.   for (i=0;  i<N2; i++) SWB(flags,swbseed[i+N-N2],swbseed[i],swbseed[i]);
  92.   for (i=N2; i<N;  i++) SWB(flags,swbseed[i  -N2],swbseed[i],swbseed[i]);
  93.   for (i=0;  i<N;  i++) *(x++) = swbseed[i] ^ (congx = congx * 69069);
  94. }
  95.  
  96. long swb32fill(rng)
  97. RNGdata *rng;
  98. { long temp;
  99.   swb32p = swb32;
  100.   SWBfill(swb32,rng);
  101.   swb32n = N-1;
  102.   return *(swb32p++);
  103. }
  104.  
  105. long  i32bit(rng)
  106. RNGdata *rng;
  107. { return  (swb32n--) ? *(swb32p++) : swb32fill(rng); }
  108.  
  109. long  i31bit(rng)
  110. RNGdata *rng;
  111. { return ((swb32n--) ? *(swb32p++) : swb32fill(rng)) & 0x7FFFFFFF; }
  112.  
  113. #define two2neg31  ( (2.0/0x10000) / 0x10000 )
  114. #define two2neg32  ( (1.0/0x10000) / 0x10000 )
  115.  
  116. /*********************************/
  117. /* External interface to mrandom */
  118. /*********************************/
  119. double _uni(rng)
  120. RNGdata *rng;
  121. {
  122.   long temp;
  123.  
  124.   temp = i31bit(rng);
  125.   if (temp & 0xFF000000) { return temp * two2neg31; }
  126.   return (temp + i32bit(rng) * two2neg32) * two2neg31;
  127. }
  128.  
  129. /* rinit initializes the constants and fills the seed array one bit at
  130.    a time by taking the leading bit of the xor of a shift register
  131.    and a congruential sequence. The same congruential generator continues
  132.    to be used as a mixing generator for the Subtract-with-borrow generator
  133.    to produce the `ultra' random numbers
  134.  
  135.    Since this is called just once, speed doesn't matter much and it might
  136.    be fine to leave this subroutine coded just as it is.
  137.  
  138.    PS:    there are quick and easy ways to fill this, but since random number
  139.     generators are really "randomness amplifiers", it is important to
  140.     start off on the right foot. This is why we take such care here.
  141. */
  142. void _rinit(rng, seed)
  143. RNGdata *rng;
  144. unsigned long *seed;
  145. {
  146.   short i,j;
  147.   unsigned long tidbits;
  148.   unsigned long congy, shrgx;
  149.  
  150.   congy=seed[0];
  151.   shrgx=seed[1];
  152.  
  153.   congx=congy*2+1;
  154.   for (i=0;i<N;i++) {
  155.     for (j=32;j>0;j--) {
  156.       congx = congx * 69069;
  157.       shrgx = shrgx ^ (shrgx >> 15);
  158.       shrgx = shrgx ^ (shrgx << 17);
  159.       tidbits = (tidbits>>1) | (0x80000000 & (congx^shrgx));
  160.     }
  161.     swbseed[i] = tidbits;
  162.   }
  163.   swb32n = 0;
  164.   flags = 0;
  165. }
  166.  
  167. int _ultra_check(rng)
  168. RNGdata *rng;
  169. {
  170. return(1);
  171. }
  172.