home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / cmath / drand.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  3.0 KB  |  156 lines
  1. /*                            drand.c
  2.  *
  3.  *    Pseudorandom number generator
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double y, drand();
  10.  *
  11.  * drand( &y );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Yields a random number 1.0 <= y < 2.0.
  18.  *
  19.  * The three-generator congruential algorithm by Brian
  20.  * Wichmann and David Hill (BYTE magazine, March, 1987,
  21.  * pp 127-8) is used. The period, given by them, is
  22.  * 6953607871644.
  23.  *
  24.  * Versions invoked by the different arithmetic compile
  25.  * time options DEC, IBMPC, and MIEEE, produce
  26.  * approximately the same sequences, differing only in the
  27.  * least significant bits of the numbers. The UNK option
  28.  * implements the algorithm as recommended in the BYTE
  29.  * article.  It may be used on all computers. However,
  30.  * the low order bits of a double precision number may
  31.  * not be adequately random, and may vary due to arithmetic
  32.  * implementation details on different computers.
  33.  *
  34.  * The other compile options generate an additional random
  35.  * integer that overwrites the low order bits of the double
  36.  * precision number.  This reduces the period by a factor of
  37.  * two but tends to overcome the problems mentioned.
  38.  *
  39.  */
  40.  
  41.  
  42.  
  43. #include "mconf.h"
  44.  
  45.  
  46. /*  Three-generator random number algorithm
  47.  * of Brian Wichmann and David Hill
  48.  * BYTE magazine, March, 1987 pp 127-8
  49.  *
  50.  * The period, given by them, is (p-1)(q-1)(r-1)/4 = 6.95e12.
  51.  * This program was written by Steve Moshier.
  52.  */
  53.  
  54. static int sx = 1;
  55. static int sy = 10000;
  56. static int sz = 3000;
  57. static double unkans = 0.0;
  58.  
  59. /* This function implements the three
  60.  * congruential generators.
  61.  */
  62.  
  63. ranwh()
  64. {
  65. int r, s;
  66.  
  67. /*  sx = sx * 171 mod 30269 */
  68. r = sx/177;
  69. s = sx - 177 * r;
  70. sx = 171 * s - 2 * r;
  71. if( sx < 0 )
  72.     sx += 30269;
  73.  
  74.  
  75. /* sy = sy * 172 mod 30307 */
  76. r = sy/176;
  77. s = sy - 176 * r;
  78. sy = 172 * s - 35 * r;
  79. if( sy < 0 )
  80.     sy += 30307;
  81.  
  82. /* sz = 170 * sz mod 30323 */
  83. r = sz/178;
  84. s = sz - 178 * r;
  85. sz = 170 * s - 63 * r;
  86. if( sz < 0 )
  87.     sz += 30323;
  88. /* The results are in static sx, sy, sz. */
  89. }
  90.  
  91. /*    drand.c
  92.  *
  93.  * Random double precision floating point number between 1 and 2.
  94.  *
  95.  * C callable:
  96.  *    drand( &x );
  97.  */
  98.  
  99. drand( a )
  100. double *a;
  101. {
  102. unsigned short r;
  103. unsigned short *p;
  104. #ifdef DEC
  105. unsigned short s, t;
  106. #endif
  107.  
  108. /* This algorithm of Wichmann and Hill computes a floating point
  109.  * result:
  110.  */
  111. ranwh();
  112. unkans = sx/30269.0  +  sy/30307.0  +  sz/30323.0;
  113. r = unkans;
  114. unkans -= r;
  115. unkans += 1.0;
  116.  
  117. /* if UNK option, do nothing further.
  118.  * Otherwise, make a random 16 bit integer
  119.  * to overwrite the least significant word
  120.  * of unkans.
  121.  */
  122. #ifdef UNK
  123. /* do nothing */
  124. #else
  125. ranwh();
  126. r = sx * sy + sz;
  127. p = (unsigned short *)&unkans;
  128. #endif
  129.  
  130. #ifdef DEC
  131. /* To make the numbers as similar as possible
  132.  * in all arithmetics, the random integer has
  133.  * to be inserted 3 bits higher up in a DEC number.
  134.  * An alternative would be put it 3 bits lower down
  135.  * in all the other number types.
  136.  */
  137. s = *(p + 2);
  138. t = s & 07;    /* save these bits to put in at the bottom */
  139. s &= 0177770;
  140. s |= (r >> 13) & 07;
  141. *(p + 2) = s;
  142. t |= r << 3;
  143. *(p + 3) = t;
  144. #endif
  145.  
  146. #ifdef IBMPC
  147. *p = r;
  148. #endif
  149.  
  150. #ifdef MIEEE
  151. *(p + 3) = r;
  152. #endif
  153.  
  154. *a = unkans;
  155. }
  156.