home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / diverses / jpl_c / random.c < prev    next >
Encoding:
C/C++ Source or Header  |  1986-05-18  |  5.1 KB  |  138 lines
  1. /* 1.1  10-16-84 */
  2. /************************************************************************
  3.  *            Robert C. Tausworthe                *
  4.  *            Jet Propulsion Laboratory            *
  5.  *            Pasadena, CA 91009        1984        *
  6.  ************************************************************************
  7.  *    Uniform Random Number Generator, programmed using the algorithm
  8.  *    given in:
  9.  *
  10.  *    Tausworthe, Robert C., STANDARDIZED DEVELOPMENT OF COMPUTER
  11.  *    SOFTWARE, Part 2, Appendix L, Example 3, Prentice-Hall, Inc.,
  12.  *    Englewood Cliffs, NJ, 1979, pp. 499-506.
  13.  *
  14.  *    The generator uses the primitive polynomial over GF(2^521),
  15.  *    f(x) = x^521 + x^32 + 1, found in Zierler, N., and Brillhart, J.,
  16.  *    "On Primitive Polynomials (Mod 2), II," INFORMATION AND CONTROL,
  17.  *    Volume 14, No. 6, pp. 566-569, June, 1969.
  18.  *
  19.  *    Numbers output by the generator are unsigned, uniformly distributed
  20.  *    over 0..65535.  Numbers have no serial statistical correlation,
  21.  *    adjacent numbers are uniformly distributed over the 32-cube, and
  22.  *    runs-up/down of length up to 13 require more than 10^11 samples to
  23.  *    show deviation from theoretical.  The period of the generator is
  24.  *    slightly larger than 1.3 * 10^154.
  25.  *
  26.  *    For more information on Tausworthe random number generators, see;
  27.  *
  28.  *    Tausworthe, Robert C., "Random Numbers Generated by Linear Recurrence
  29.  *    Modulo 2," MATH COMP., Vol. XIX, No. 90, pp. 201-208, April 1965.
  30.  *
  31.  *    Toothill, J. P. R., "An Asymptotically Random Tausworthe Sequence,"
  32.  *    ACM J., Vol. 20, No. 3, pp. 469-481, July 1973.
  33.  *
  34.  *    Toothill, J. P. R., et al., "The Runs-Up and -Down Performance of
  35.  *    Tausworthe Pseudorandom Number Generators," ACM J., Vol. 18, No. 3,
  36.  *    pp. 381-399, July 1971.
  37.  *
  38.  *----------------------------------------------------------------------*
  39.  *
  40.  *    The shift register of 521 bits requires 65.125 bytes, or 65 bytes
  41.  *    plus one bit.  One register load then provides 32 16-bit numbers
  42.  *    before it must be refilled.  The 32-bit tap is at a 4-byte shift
  43.  *    from the first byte, and at 521 - 32 = 489 = (61 bytes + 1 bit)
  44.  *    inverse shift from the end.
  45.  *
  46.  *----------------------------------------------------------------------*/
  47.  
  48. #include "defs.h"
  49. #include "stdtyp.h"
  50.  
  51. /*----------------------------------------------------------------------*/
  52.  
  53. #define MULTIPLIER    18829    /* This is 5^15 (mod 65536)        */
  54. #define PRESEED        0xe5a1    /* a number picked out of the air!    */
  55.  
  56. /*\p--------------------------------------------------------------------*/
  57.  
  58. LOCAL utiny shiftreg[66];        /* 521-bit shift register    */
  59. LOCAL utiny *randptr = shiftreg - 1;    /* initialized for not accessed    */
  60. LOCAL utiny *endreg = shiftreg + 62;    /* always points at 32nd number    */
  61.  
  62. /************************************************************************/
  63.     VOID
  64. randize(seed)            /* Randomize the generator using seed.
  65.                    See algorithm for seeding values.    */
  66. /*----------------------------------------------------------------------*/
  67. {
  68.     int i, *rseed;
  69.  
  70.     if (seed IS 0)    
  71.         seed = PRESEED;
  72.     else if (seed < 0)
  73.     {    rseed = (int *) &shiftreg[16];
  74.         while ((seed = *(++rseed)) IS 0)
  75.             ;    /* find something non-zero, whatever.    */
  76.     }
  77.     for (i = 0; i < 66; i++)
  78.         /* assume no adverse reaction on integer overflow    */
  79.     {    shiftreg[i++] = ((seed *= MULTIPLIER) &0xff);
  80.         shiftreg[i] = ((seed >> 8) & 0xff);
  81.     }
  82.     shiftreg[65] &= 0x80;    /* mask out bits beyond 521st.        */
  83.     refill();    /* refill and reset randptr            */
  84. }
  85.  
  86. /************************************************************************/
  87.     unsigned
  88. urandom()    /* Return the next random number.  If shift register is
  89.            depleted, refill the entire register.        */
  90. /*----------------------------------------------------------------------*/
  91. {
  92.     unsigned r;
  93.  
  94.     if (randptr < shiftreg)
  95.         randize(0);    /* in case randize wasn't called first    */
  96.     else if (randptr > endreg)
  97.         refill();    /* and reset randptr & randx        */
  98.     r = UTINY(*randptr++);
  99.     r += (UTINY(*randptr++) << 8);
  100.     return (r);
  101. }
  102. /*\p*********************************************************************/
  103.     double
  104. random()        /* Return a uniform random value in [0, 1].    */
  105.  
  106. /*----------------------------------------------------------------------*/
  107. {
  108.     unsigned urandom();
  109.  
  110.     return (urandom() / 65535.0);
  111. }
  112.  
  113. /************************************************************************/
  114.     LOCAL VOID
  115. refill()            /* refill the shift register        */
  116.     
  117. /*----------------------------------------------------------------------*/
  118. {
  119.     utiny *p, *q, cy0, cy1;
  120.     FAST int i;
  121.  
  122.     p = shiftreg;            /* point at 1st byte.        */
  123.     q = &p[4];            /* shift 4 * 8 = 32,        */
  124.     for (i = 0; i < 62; i++)    /* and mod-2 add the register    */
  125.         *p++ ^= *q++;        /* to the 32-bit shift.        */
  126.     p--;                /* p should now be at byte 61:    */
  127.     q = shiftreg;            /* first 489 bits are ok.    */
  128.     cy0 = 0;            /* set carry bit to zero    */
  129.     for (i = 0; i < 4; i++)        /* for the remaining 32 bits:    */
  130.     {    cy1 = ((*q & 1) ? 0x80 : 0); /* save carry for nxt shift*/
  131.         *p++ ^= (((*q++ >> 1) & 0x7f) | cy0);
  132.                     /* mod-2 add the 489-shifts    */
  133.         cy0 = cy1;        /* set carry for next shift    */
  134.     }
  135.     *p ^= cy0;            /* and mod-2 the final bit.    */
  136.     randptr = shiftreg;        /* point at the first number    */
  137. }
  138.