home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / cprob / drand.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  2.9 KB  |  155 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.  */
  52.  
  53. static int sx = 1;
  54. static int sy = 10000;
  55. static int sz = 3000;
  56. static double unkans = 0.0;
  57.  
  58. /* This function implements the three
  59.  * congruential generators.
  60.  */
  61.  
  62. ranwh()
  63. {
  64. int r, s;
  65.  
  66. /*  sx = sx * 171 mod 30269 */
  67. r = sx/177;
  68. s = sx - 177 * r;
  69. sx = 171 * s - 2 * r;
  70. if( sx < 0 )
  71.     sx += 30269;
  72.  
  73.  
  74. /* sy = sy * 172 mod 30307 */
  75. r = sy/176;
  76. s = sy - 176 * r;
  77. sy = 172 * s - 35 * r;
  78. if( sy < 0 )
  79.     sy += 30307;
  80.  
  81. /* sz = 170 * sz mod 30323 */
  82. r = sz/178;
  83. s = sz - 178 * r;
  84. sz = 170 * s - 63 * r;
  85. if( sz < 0 )
  86.     sz += 30323;
  87. /* The results are in static sx, sy, sz. */
  88. }
  89.  
  90. /*    drand.c
  91.  *
  92.  * Random double precision floating point number between 1 and 2.
  93.  *
  94.  * C callable:
  95.  *    drand( &x );
  96.  */
  97.  
  98. drand( a )
  99. double *a;
  100. {
  101. unsigned short r;
  102. unsigned short *p;
  103. #ifdef DEC
  104. unsigned short s, t;
  105. #endif
  106.  
  107. /* This algorithm of Wichmann and Hill computes a floating point
  108.  * result:
  109.  */
  110. ranwh();
  111. unkans = sx/30269.0  +  sy/30307.0  +  sz/30323.0;
  112. r = unkans;
  113. unkans -= r;
  114. unkans += 1.0;
  115.  
  116. /* if UNK option, do nothing further.
  117.  * Otherwise, make a random 16 bit integer
  118.  * to overwrite the least significant word
  119.  * of unkans.
  120.  */
  121. #ifdef UNK
  122. /* do nothing */
  123. #else
  124. ranwh();
  125. r = sx * sy + sz;
  126. p = (unsigned short *)&unkans;
  127. #endif
  128.  
  129. #ifdef DEC
  130. /* To make the numbers as similar as possible
  131.  * in all arithmetics, the random integer has
  132.  * to be inserted 3 bits higher up in a DEC number.
  133.  * An alternative would be put it 3 bits lower down
  134.  * in all the other number types.
  135.  */
  136. s = *(p + 2);
  137. t = s & 07;    /* save these bits to put in at the bottom */
  138. s &= 0177770;
  139. s |= (r >> 13) & 07;
  140. *(p + 2) = s;
  141. t |= r << 3;
  142. *(p + 3) = t;
  143. #endif
  144.  
  145. #ifdef IBMPC
  146. *p = r;
  147. #endif
  148.  
  149. #ifdef MIEEE
  150. *(p + 3) = r;
  151. #endif
  152.  
  153. *a = unkans;
  154. }
  155.