home *** CD-ROM | disk | FTP | other *** search
/ Serving the Web / ServingTheWeb1995.disc1of1.iso / linux / slacksrce / d / libc / libc-4.6 / libc-4 / libc-linux / misc / drand48.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-02-28  |  4.1 KB  |  185 lines
  1. /*    @(#)drand48.c    2.2    */
  2. /*LINTLIBRARY*/
  3. /*
  4.  *    drand48, etc. pseudo-random number generator
  5.  *    This implementation assumes unsigned short integers of at least
  6.  *    16 bits, long integers of at least 32 bits, and ignores
  7.  *    overflows on adding or multiplying two unsigned integers.
  8.  *    Two's-complement representation is assumed in a few places.
  9.  *    Some extra masking is done if unsigneds are exactly 16 bits
  10.  *    or longs are exactly 32 bits, but so what?
  11.  *    An assembly-language implementation would run significantly faster.
  12.  */
  13.  
  14. #include <stdlib.h>
  15.  
  16. #ifndef HAVEFP
  17. #define HAVEFP 1
  18. #endif
  19. #define N    16
  20. #define MASK    ((unsigned)(1 << (N - 1)) + (1 << (N - 1)) - 1)
  21. #define LOW(x)    ((unsigned)(x) & MASK)
  22. #define HIGH(x)    LOW((x) >> N)
  23. #define MUL(x, y, z)    { long l = (long)(x) * (long)(y); \
  24.         (z)[0] = LOW(l); (z)[1] = HIGH(l); }
  25. #define CARRY(x, y)    ((long)(x) + (long)(y) > MASK)
  26. #define ADDEQU(x, y, z)    (z = CARRY(x, (y)), x = LOW(x + (y)))
  27. #define X0    0x330E
  28. #define X1    0xABCD
  29. #define X2    0x1234
  30. #define A0    0xE66D
  31. #define A1    0xDEEC
  32. #define A2    0x5
  33. #define C    0xB
  34. #define SET3(x, x0, x1, x2)    ((x)[0] = (x0), (x)[1] = (x1), (x)[2] = (x2))
  35. #define SETLOW(x, y, n) SET3(x, LOW((y)[n]), LOW((y)[(n)+1]), LOW((y)[(n)+2]))
  36. #define SEED(x0, x1, x2) (SET3(x, x0, x1, x2), SET3(a, A0, A1, A2), c = C)
  37. #define REST(v)    for (i = 0; i < 3; i++) { xsubi[i] = x[i]; x[i] = temp[i]; } \
  38.         return (v);
  39. #define NEST(TYPE, f, F)    TYPE f(xsubi) register unsigned short int *xsubi; { \
  40.     register int i; register TYPE v; unsigned temp[3]; \
  41.     for (i = 0; i < 3; i++) { temp[i] = x[i]; x[i] = LOW(xsubi[i]); }  \
  42.     v = F(); REST(v); }
  43. #define HI_BIT    (1L << (2 * N - 1))
  44.  
  45. static unsigned x[3] = { X0, X1, X2 }, a[3] = { A0, A1, A2 }, c = C;
  46. static unsigned short lastx[3];
  47. static void next();
  48.  
  49. #if HAVEFP
  50. double
  51. drand48()
  52. {
  53. #if pdp11
  54.     static double two16m; /* old pdp11 cc can't compile an expression */
  55.     two16m = 1.0 / (1L << N); /* in "double" initializer! */
  56. #else
  57.     static double two16m = 1.0 / (1L << N);
  58. #endif
  59.  
  60.     next();
  61.     return (two16m * (two16m * (two16m * x[0] + x[1]) + x[2]));
  62. }
  63.  
  64. NEST(double, erand48, drand48);
  65.  
  66. #else
  67.  
  68. long
  69. irand48(m)
  70. /* Treat x[i] as a 48-bit fraction, and multiply it by the 16-bit
  71.  * multiplier m.  Return integer part as result.
  72.  */
  73. register unsigned short m;
  74. {
  75.     unsigned r[4], p[2], carry0 = 0;
  76.  
  77.     next();
  78.     MUL(m, x[0], &r[0]);
  79.     MUL(m, x[2], &r[2]);
  80.     MUL(m, x[1], p);
  81.     if (CARRY(r[1], p[0]))
  82.         ADDEQU(r[2], 1, carry0);
  83.     return (r[3] + carry0 + CARRY(r[2], p[1]));
  84. }
  85.  
  86. long
  87. krand48(xsubi, m)
  88. /* same as irand48, except user provides storage in xsubi[] */
  89. register unsigned short *xsubi;
  90. unsigned short m;
  91. {
  92.     register int i;
  93.     register long iv;
  94.     unsigned temp[3];
  95.  
  96.     for (i = 0; i < 3; i++) {
  97.         temp[i] = x[i];
  98.         x[i] = xsubi[i];
  99.     }
  100.     iv = irand48(m);
  101.     REST(iv);
  102. }
  103. #endif
  104.  
  105. long int
  106. lrand48()
  107. {
  108.     next();
  109.     return (((long)x[2] << (N - 1)) + (x[1] >> 1));
  110. }
  111.  
  112. long int
  113. mrand48()
  114. {
  115.     register long l;
  116.  
  117.     next();
  118.     /* sign-extend in case length of a long > 32 bits
  119.                         (as on Honeywell) */
  120.     return ((l = ((long)x[2] << N) + x[1]) & HI_BIT ? l | -HI_BIT : l);
  121. }
  122.  
  123. static void
  124. next()
  125. {
  126.     unsigned p[2], q[2], r[2], carry0, carry1;
  127.  
  128.     MUL(a[0], x[0], p);
  129.     ADDEQU(p[0], c, carry0);
  130.     ADDEQU(p[1], carry0, carry1);
  131.     MUL(a[0], x[1], q);
  132.     ADDEQU(p[1], q[0], carry0);
  133.     MUL(a[1], x[0], r);
  134.     x[2] = LOW(carry0 + carry1 + CARRY(p[1], r[0]) + q[1] + r[1] +
  135.         a[0] * x[2] + a[1] * x[1] + a[2] * x[0]);
  136.     x[1] = LOW(p[1] + r[0]);
  137.     x[0] = LOW(p[0]);
  138. }
  139.  
  140. void
  141. srand48(long int seedval)
  142. {
  143.     SEED(X0, LOW(seedval), HIGH(seedval));
  144. }
  145.  
  146. unsigned short int *
  147. seed48(unsigned short seed16v[3])
  148. {
  149.     SETLOW(lastx, x, 0);
  150.     SEED(LOW(seed16v[0]), LOW(seed16v[1]), LOW(seed16v[2]));
  151.     return (lastx);
  152. }
  153.  
  154. void
  155. lcong48(unsigned short int param[7])
  156. {
  157.     SETLOW(x, param, 0);
  158.     SETLOW(a, param, 3);
  159.     c = LOW(param[6]);
  160. }
  161.  
  162. NEST(long, nrand48, lrand48);
  163.  
  164. NEST(long, jrand48, mrand48);
  165.  
  166. #ifdef DRIVER
  167. /*
  168.     This should print the sequences of integers in Tables 2
  169.         and 1 of the TM:
  170.     1623, 3442, 1447, 1829, 1305, ...
  171.     657EB7255101, D72A0C966378, 5A743C062A23, ...
  172.  */
  173. #include <stdio.h>
  174.  
  175. main()
  176. {
  177.     int i;
  178.  
  179.     for (i = 0; i < 80; i++) {
  180.         printf("%4d ", (int)(4096 * drand48()));
  181.         printf("%.4X%.4X%.4X\n", x[2], x[1], x[0]);
  182.     }
  183. }
  184. #endif
  185.