home *** CD-ROM | disk | FTP | other *** search
/ Super PC 34 / Super PC 34 (Shareware).iso / spc / UTIL / DJGPP2 / V2 / DJTST200.ZIP / tests / libc / ansi / math / cephes / drand.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-15  |  3.4 KB  |  186 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. void ranwh(void);
  60.  
  61. /* This function implements the three
  62.  * congruential generators.
  63.  */
  64.  
  65. void
  66.  
  67. ranwh()
  68. {
  69. int r, s;
  70.  
  71. /*  sx = sx * 171 mod 30269 */
  72. r = sx/177;
  73. s = sx - 177 * r;
  74. sx = 171 * s - 2 * r;
  75. if( sx < 0 )
  76.     sx += 30269;
  77.  
  78.  
  79. /* sy = sy * 172 mod 30307 */
  80. r = sy/176;
  81. s = sy - 176 * r;
  82. sy = 172 * s - 35 * r;
  83. if( sy < 0 )
  84.     sy += 30307;
  85.  
  86. /* sz = 170 * sz mod 30323 */
  87. r = sz/178;
  88. s = sz - 178 * r;
  89. sz = 170 * s - 63 * r;
  90. if( sz < 0 )
  91.     sz += 30323;
  92. /* The results are in static sx, sy, sz. */
  93. }
  94.  
  95. /*    drand.c
  96.  *
  97.  * Random double precision floating point number between 1 and 2.
  98.  *
  99.  * C callable:
  100.  *    drand( &x );
  101.  */
  102.  
  103. void
  104. drandl(long double *a)
  105. {
  106.     double x;
  107.     drand(&x);
  108.     *a = (long double)x;
  109. }
  110.  
  111. void
  112. drand(double  *a)
  113. {
  114. unsigned short r;
  115. unsigned short *p;
  116. #ifdef DEC
  117. unsigned short s, t;
  118. #endif
  119.  
  120. /* This algorithm of Wichmann and Hill computes a floating point
  121.  * result:
  122.  */
  123. ranwh();
  124. unkans = sx/30269.0  +  sy/30307.0  +  sz/30323.0;
  125. r = unkans;
  126. unkans -= r;
  127. unkans += 1.0;
  128.  
  129. /* if UNK option, do nothing further.
  130.  * Otherwise, make a random 16 bit integer
  131.  * to overwrite the least significant word
  132.  * of unkans.
  133.  */
  134. #ifdef UNK
  135. /* do nothing */
  136. #else
  137. ranwh();
  138. r = sx * sy + sz;
  139. p = (unsigned short *)&unkans;
  140. #endif
  141.  
  142. #ifdef DEC
  143. /* To make the numbers as similar as possible
  144.  * in all arithmetics, the random integer has
  145.  * to be inserted 3 bits higher up in a DEC number.
  146.  * An alternative would be put it 3 bits lower down
  147.  * in all the other number types.
  148.  */
  149. s = *(p + 2);
  150. t = s & 07;    /* save these bits to put in at the bottom */
  151. s &= 0177770;
  152. s |= (r >> 13) & 07;
  153. *(p + 2) = s;
  154. t |= r << 3;
  155. *(p + 3) = t;
  156. #endif
  157.  
  158. #ifdef IBMPC
  159. *p = r;
  160. #endif
  161.  
  162. #ifdef MIEEE
  163. *(p + 3) = r;
  164. #endif
  165.  
  166. *a = unkans;
  167. }
  168.  
  169. #ifdef TEST
  170. #include <stdio.h>
  171.  
  172.  
  173. int
  174. main()
  175. {
  176.     double d;
  177.     int i;
  178.     for (i=0; i<100; i++)
  179.     {
  180.         drand(&d);
  181.         printf("%.20g\n", d);
  182.     }
  183.     return(0);
  184. }
  185. #endif
  186.