home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / ixemul-45.0-src.tgz / tar.out / contrib / ixemul / static / rand48.c < prev    next >
C/C++ Source or Header  |  1996-10-01  |  6KB  |  233 lines
  1. /* Random number generation using the linear congruential algorithm
  2.    X(n+1) = (a * X(n) + c) mod m
  3.    with a precision of 48 bits.
  4.  
  5.    Author: Kriton Kyrimis (kyrimis@cti.gr)
  6.    Code status: Public Domain.
  7. */
  8.  
  9. #include <limits.h>
  10. #include <float.h>
  11.  
  12. /* Parameters for the linear congruential algorithm:
  13.    parm[0..2] is the current value of Xn (internal seed, m.s.word last)
  14.    parm[3..5] is the value of a (m.s.word last)
  15.    parm[6] is the value of c.
  16. */
  17. #define X0 0x1234 /* MSB * Initial value for Xn, obtained using seed48() */
  18. #define X1 0xABCD    /* on SunOS 4.1.3 */
  19. #define X2 0x330E
  20.  
  21. #define A0 0x0005 /* MSB * Default value for a, taken from the man page */
  22. #define A1 0xDEEC
  23. #define A2 0xE66D
  24.  
  25. #define C0 0x000B    /* Default value for c, taken from the man page */
  26.  
  27. static unsigned short parm[7] = {
  28.   X2, X1, X0,
  29.   A2, A1, A0,
  30.   C0
  31. };
  32.  
  33. /* To produce a double random number in [0,1) we get a 32-bit unsigned long
  34.    random number, convert it to double, and divide it by ULONG_MAX + EPSILON.
  35.    (We add the EPSILON to exclude 1.0 from the set of possible results.) 
  36.  
  37.    We derive EPSILON by noting that for a random value of ULONG_MAX,
  38.    we want to return the smallest double that is less than 1.0.
  39.    Therefore:
  40.  
  41.         ULONG_MAX
  42.    --------------------- = (1.0 - DBL_EPSILON)
  43.    (ULONG_MAX + EPSILON)
  44.  
  45.    (This is probably overkill.)
  46.  
  47. */
  48.  
  49. #define EPSILON    (double)ULONG_MAX*(1.0/(1.0-DBL_EPSILON)-1.0)
  50.  
  51.  
  52. /*--------------------------------------------------------------------------*
  53.  * Parameter initialization functions                                       *
  54.  *--------------------------------------------------------------------------*/
  55.  
  56. /* This function sets the two m.s.words of the internal seed to the value
  57.    supplied by the caller, and the l.s.word of the internal seed to 0x330E.
  58. */
  59. void
  60. srand48(long seed)
  61. {
  62.   parm[0] = 0x330E;
  63.   parm[1] = ((unsigned long)seed) & 0xFFFF;
  64.   parm[2] = ((unsigned long)seed >> 16) & 0xFFFF;
  65.   parm[3] = A2;
  66.   parm[4] = A1;
  67.   parm[5] = A0;
  68.   parm[6] = C0;
  69. }
  70.  
  71. /* This function sets all three words of the internal seed to the value
  72.    supplied by the caller. It returns a pointer to an array containing
  73.    a copy of the previous value of the internal seed.
  74. */
  75. unsigned short *
  76. seed48(unsigned short *seed)
  77. {
  78.   static unsigned short oldparm[3];
  79.   unsigned short tmpparm[3];
  80.  
  81.   /* Can't assign oldparm[] = parm[] directly, because seed[] may be a pointer
  82.      to oldparm[], obtained from a previous call to seed48 , in which case
  83.      we would destroy the contents of seed[] */
  84.   tmpparm[0] = parm[0];
  85.   tmpparm[1] = parm[1];
  86.   tmpparm[2] = parm[2];
  87.   parm[0] = seed[0];
  88.   parm[1] = seed[1];
  89.   parm[2] = seed[2];
  90.   oldparm[0] = tmpparm[0];
  91.   oldparm[1] = tmpparm[1];
  92.   oldparm[2] = tmpparm[2];
  93.   parm[3] = A2;
  94.   parm[4] = A1;
  95.   parm[5] = A0;
  96.   parm[6] = C0;
  97.  
  98.   return oldparm;
  99. }
  100.  
  101. /* This function sets all seven words of the internal parameters array to the
  102.    values supplied by the caller.
  103. */
  104. void
  105. lcong48(unsigned short *new_parm)
  106. {
  107.   parm[0] = new_parm[0];
  108.   parm[1] = new_parm[1];
  109.   parm[2] = new_parm[2];
  110.   parm[3] = new_parm[3];
  111.   parm[4] = new_parm[4];
  112.   parm[5] = new_parm[5];
  113.   parm[6] = new_parm[6];
  114. }
  115.  
  116.  
  117. /*--------------------------------------------------------------------------*
  118.  * Random number generator                                                  *
  119.  *--------------------------------------------------------------------------*/
  120.  
  121. /* This function implements the linear congruential algorithm.  Thanks to
  122.    gcc's long long ints, implementing 48-bit arithmetic (actually 64-bit,
  123.    truncating the result) is trivial.  Limitations of long long int
  124.    implementation in (amiga?) gcc 2.7.0, made me use the kludge with the union
  125.    to convert from short[3] to long long int. (It's probably faster though!)
  126.  
  127.    This function takes an array of three shorts (a 48-bit seed, m.s.word
  128.    last) and returns a long between -2**31 and 2**31-1, updating the seed
  129.    (the result is the two m.s.words of the updated seed).
  130. */
  131.  
  132. static long
  133. rng(unsigned short *seed)
  134. {
  135.   long long int Xn, Xn1, a, c;
  136.   union {
  137.     long long int l;
  138.     unsigned short s[4];
  139.   } i;
  140.  
  141.   i.s[0] = 0;
  142.   i.s[1] = seed[2];
  143.   i.s[2] = seed[1];
  144.   i.s[3] = seed[0];
  145.   Xn = i.l;
  146.  
  147.   i.s[0] = 0;
  148.   i.s[1] = parm[5];
  149.   i.s[2] = parm[4];
  150.   i.s[3] = parm[3];
  151.   a = i.l;
  152.  
  153.   c =  (long long int)(parm[6]);
  154.  
  155.   Xn1 = a * Xn + c;
  156.  
  157.   i.l = Xn1;
  158.   seed[0] = i.s[3];
  159.   seed[1] = i.s[2];
  160.   seed[2] = i.s[1];
  161.  
  162.   return (long)((((unsigned long)seed[2]) << 16) + seed[1]);
  163. }
  164.  
  165.  
  166. /*--------------------------------------------------------------------------*
  167.  * Interface functions to the random number generator                       *
  168.  *--------------------------------------------------------------------------*/
  169.  
  170. /* This function returns a long between 0 and 2**31-1 by calling rng
  171.    with the internal seed, returning the 31 most significant bits of the
  172.    result shifted by one position to the right.
  173. */
  174. long
  175. lrand48(void)
  176. {
  177.   return (rng(parm) >> 1) & 0x7FFFFFFF;
  178. }
  179.  
  180. /* Same as lrand48(), but using an external seed. */
  181.  
  182. long
  183. nrand48(unsigned short seed[3])
  184. {
  185.   return (rng(seed) >> 1) & 0x7FFFFFFF;
  186. }
  187.  
  188. /* This function returns a long between -2**31 and 2**31-1 by calling rng
  189.    with the internal seed.
  190. */
  191. long
  192. mrand48(void)
  193. {
  194.   return rng(parm);
  195. }
  196.  
  197. /* Same as mrand48(), but using an external seed. */
  198.  
  199. long
  200. jrand48(unsigned short seed[3])
  201. {
  202.   return rng(seed);
  203. }
  204.  
  205. /* This function returns a double in the interval [0,1) by calling mrand48()
  206.    and dividing the result by ULONG_MAX + EPSILON. */
  207.  
  208. double
  209. drand48(void)
  210. {
  211.   union {
  212.     long l;
  213.     unsigned long u;
  214.   } x;
  215.  
  216.   x.l = mrand48();
  217.   return (double)x.u / ((double)(ULONG_MAX) + EPSILON);
  218. }
  219.  
  220. /* Same as drand48(), but using an external seed. */
  221.  
  222. double
  223. erand48(unsigned short seed[3])
  224. {
  225.   union {
  226.     long l;
  227.     unsigned long u;
  228.   } x;
  229.  
  230.   x.l = nrand48(seed);
  231.   return (double)x.u / ((double)(ULONG_MAX) + EPSILON);
  232. }
  233.