home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume28 / ldb / part02 / r_xrand.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-03-15  |  8.4 KB  |  269 lines
  1. /*    r_xrand.c
  2.  */
  3.  
  4. #include "ldb.h"
  5. #include <math.h>
  6.  
  7. /*======================================================================
  8.  * This file contains the random number generator.  It is accessed by
  9.  * calling RandomInit with a 32-bit seed, then calling Rolldie() to
  10.  * generate integers in the range [1-6].  This particular implementation
  11.  * was taken from xrand in volume 2 of comp.sources.misc.  It was
  12.  * written by Andreas Nowatzyk, then of Carnegie-Mellon University, and
  13.  * is reproduced here essentially unchanged (although I omitted
  14.  * all routines but the ones needed for integers in [1-6]).
  15.  * It is used by permission of the author, whose original copyright
  16.  * and permission-to-use appears below.
  17.  *======================================================================
  18.  */
  19.  
  20.  
  21.  
  22. /*----------------------------------------------------------------------
  23.  *    Rolldie -- return a random number between 1 and 6
  24.  *----------------------------------------------------------------------
  25.  */
  26.  
  27. Rolldie()
  28. {
  29.  
  30. return(rnd_ri(6)+1);
  31. }
  32.  
  33.  
  34.  
  35. /*----------------------------------------------------------------------
  36.  *    RandomInit -- initialize the random number generator
  37.  *
  38.  * This function should be called once, before Rolldie is called.
  39.  * The seed may be any more or less random number; the time seems
  40.  * to be conventional.
  41.  *----------------------------------------------------------------------
  42.  */
  43.  
  44. RandomInit(seed)
  45. long seed;
  46. {
  47.  
  48. rnd_init(seed);
  49. }
  50.  
  51.  
  52.  
  53. /*----------------------------------------------------------------------
  54.  *From:    EMF780::WINS%"Andreas.Nowatzyk@eng.sun.com" 20-NOV-1991 17:50:38.73
  55.  *To:    ROSS
  56.  *CC:    
  57.  *Subj:    Re:  xrand
  58.  *
  59.  *Return-Path: <Andreas.Nowatzyk@eng.sun.com>
  60.  *Received: from everest.den.mmc.com by emf780.den.mmc.com with SMTP ; 
  61.  *          Wed, 20 Nov 91 16:50:29 MDT
  62.  *Received: from Sun.COM by everest.den.mmc.com (4.1/1.34.a)
  63.  *    id AA27426; Wed, 20 Nov 91 17:52:10 MST
  64.  *Received: from Eng.Sun.COM (zigzag-bb.Corp.Sun.COM) by Sun.COM (4.1/SMI-4.1)
  65.  *    id AA07190; Wed, 20 Nov 91 16:50:53 PST
  66.  *Received: from bovic.Eng.Sun.COM by Eng.Sun.COM (4.1/SMI-4.1)
  67.  *    id AA21474; Wed, 20 Nov 91 16:49:44 PST
  68.  *Received: by bovic.Eng.Sun.COM (4.1/SMI-4.1)
  69.  *    id AA02462; Wed, 20 Nov 91 16:49:16 PST
  70.  *Date: Wed, 20 Nov 91 16:49:16 PST
  71.  *From: Andreas.Nowatzyk@eng.sun.com (Andreas G. Nowatzyk)
  72.  *Message-Id: <9111210049.AA02462@bovic.Eng.Sun.COM>
  73.  *To: ross@emf780.den.mmc.com
  74.  *Subject: Re:  xrand
  75.  *
  76.  *I have no objection to distributing xrand, provided that this is not done
  77.  *for profit and that the source is credited.
  78.  *
  79.  *  --  Andreas
  80.  *
  81.  *PS: There was a bug is some version of xrand.c: The costant for the most
  82.  *    positive integer (0x7fffffff) had an 'f' missing in two places. Its should
  83.  *    have used a #define MAX_INT with apropriate value instead of hardwiring
  84.  *    an hex-constant.
  85.  *----------------------------------------------------------------------
  86.  */
  87.  
  88. #define MAX_INT 0x7fffffff
  89.  
  90. /* Random number generators:
  91.  *
  92.  *  rnd_init (unsigned seed) 
  93.  *            : initializes the generator
  94.  *
  95.  *  rnd_i ()        : returns positive integers [0,0x7fffffff]
  96.  *  rnd_ri (long n)    : returns positive integers [0,n-1]
  97.  *
  98.  *  Algorithm M as describes in Knuth's "Art of Computer Programming", Vol 2. 1969
  99.  *  is used with a linear congruential generator (to get a good uniform
  100.  *  distribution) that is permuted with a Fibonacci additive congruential
  101.  *  generator to get good independence.
  102.  *
  103.  *  Bit, byte, and word distributions were extensively tested and pass
  104.  *  Chi-squared test near perfect scores (>7E8 numbers tested, Uniformity
  105.  *  assumption holds with probability > 0.999)
  106.  *
  107.  *  Run-up tests for on 7E8 numbers confirm independence with
  108.  *  probability > 0.97.
  109.  *
  110.  *  Plotting random points in 2d reveals no apparent structure.
  111.  *
  112.  *  Autocorrelation on sequences of 5E5 numbers (A(i) = SUM X(n)*X(n-i), i=1..512)
  113.  *  results in no obvious structure (A(i) ~ const).
  114.  *
  115.  *  On a SUN 3/60, rnd_u() takes about 19.4 usec per call, which is about 44%
  116.  *  slower than Berkeley's random() (13.5 usec/call).
  117.  *
  118.  *  Except for speed and memory requirements, this generator outperforms
  119.  *  random() for all tests. (random() scored rather low on uniformity tests,
  120.  *  while independence test differences were less dramatic).
  121.  *
  122.  *  Thanks to M.Mauldin, H.Walker, J.Saxe and M.Molloy for inspiration & help.
  123.  *
  124.  *  (c) Copyright 1988 by A. Nowatzyk
  125.  *
  126.  */
  127.  
  128. /* LC-parameter selection follows recommendations in 
  129.  * "Handbook of Mathematical Functions" by Abramowitz & Stegun 10th, edi.
  130.  */
  131. #define LC_A 66049            /* = 251^2, ~= sqrt(2^32)            */
  132. #define LC_C 3907864577            /* result of a long trial & error series    */
  133.  
  134. #define Xrnd(x) (x * LC_A + LC_C)   /* the LC polynomial            */
  135.             
  136. static unsigned long Fib[55];        /* will use X(n) = X(n-55) - X(n-24)    */
  137. static int Fib_ind;            /* current index in circular buffer        */
  138. static unsigned long Xrnd_var;        /* LCA - recurrence variable        */
  139. static unsigned long auxtab[256];   /* temporal permutation table        */
  140. static unsigned long prmtab[64] = { /* spatial permutation table        */
  141.     0xffffffff, 0x00000000,  0x00000000,  0x00000000,  /* 3210 */
  142.     0x0000ffff, 0x00ff0000,  0x00000000,  0xff000000,  /* 2310 */
  143.     0xff0000ff, 0x0000ff00,  0x00000000,  0x00ff0000,  /* 3120 */
  144.     0x00ff00ff, 0x00000000,  0xff00ff00,  0x00000000,  /* 1230 */
  145.  
  146.     0xffff0000, 0x000000ff,  0x00000000,  0x0000ff00,  /* 3201 */
  147.     0x00000000, 0x00ff00ff,  0x00000000,  0xff00ff00,  /* 2301 */
  148.     0xff000000, 0x00000000,  0x000000ff,  0x00ffff00,  /* 3102 */
  149.     0x00000000, 0x00000000,  0x00000000,  0xffffffff,  /* 2103 */
  150.  
  151.     0xff00ff00, 0x00000000,  0x00ff00ff,  0x00000000,  /* 3012 */
  152.     0x0000ff00, 0x00000000,  0x00ff0000,  0xff0000ff,  /* 2013 */
  153.     0x00000000, 0x00000000,  0xffffffff,  0x00000000,  /* 1032 */
  154.     0x00000000, 0x0000ff00,  0xffff0000,  0x000000ff,  /* 1023 */
  155.  
  156.     0x00000000, 0xffffffff,  0x00000000,  0x00000000,  /* 0321 */
  157.     0x00ffff00, 0xff000000,  0x00000000,  0x000000ff,  /* 0213 */
  158.     0x00000000, 0xff000000,  0x0000ffff,  0x00ff0000,  /* 0132 */
  159.     0x00000000, 0xff00ff00,  0x00000000,  0x00ff00ff   /* 0123 */
  160. };
  161.  
  162. union hack {                /* used to access doubles as unsigneds    */
  163.     double d;
  164.     unsigned long u[2];
  165. };
  166.  
  167. static union hack man;            /* mantissa bit vector            */
  168.  
  169. rnd_init (seed)                /* modified: seed 0-31 use precomputed stuff */
  170.     unsigned seed;
  171. {
  172.     register unsigned long u;
  173.     register int i;
  174.     double x, y;
  175.     union hack t;
  176.  
  177.     static unsigned seed_tab[32] = {
  178.         0xbdcc47e5, 0x54aea45d, 0xec0df859, 0xda84637b,
  179.         0xc8c6cb4f, 0x35574b01, 0x28260b7d, 0x0d07fdbf,
  180.         0x9faaeeb0, 0x613dd169, 0x5ce2d818, 0x85b9e706,
  181.         0xab2469db, 0xda02b0dc, 0x45c60d6e, 0xffe49d10,
  182.         0x7224fea3, 0xf9684fc9, 0xfc7ee074, 0x326ce92a,
  183.         0x366d13b5, 0x17aaa731, 0xeb83a675, 0x7781cb32,
  184.         0x4ec7c92d, 0x7f187521, 0x2cf346b4, 0xad13310f,
  185.         0xb89cff2b, 0x12164de1, 0xa865168d, 0x32b56cdf  };
  186.  
  187.     if (seed < 32)
  188.     u = seed_tab[seed];
  189.     else
  190.     u = seed ^ seed_tab[seed & 31];
  191.  
  192.     for (i = 55; i--;)            /* set up Fibonacci additive congruential    */
  193.     Fib[i] = u = Xrnd(u);
  194.  
  195.     for (i = 256; i--;)
  196.     auxtab[i] = u = Xrnd(u);
  197.  
  198.     Fib_ind = u % 55;            /* select a starting point            */
  199.  
  200.     Xrnd_var = u;
  201.  
  202.     if (sizeof(x) != 2 * sizeof(unsigned long)) {
  203.     x = 0.0;
  204.     y = 1.0;
  205.     y /= x;                /*** intentional divide by 0: rnd_01d will
  206.                      not work because a double doesn't fit
  207.                      in 2 unsigned longs on your machine! ***/
  208.     };
  209.  
  210.     x = 1.0;
  211.     y = 0.5;
  212.     do {                /* find largest fp-number < 2.0        */
  213.     t.d = x;
  214.     x += y;
  215.     y *= 0.5;
  216.     } while (x != t.d && x < 2.0);
  217.  
  218.     man.d = 1.0;
  219.     man.u[0] ^= t.u[0];
  220.     man.u[1] ^= t.u[1];            /* man is now 1 for each mantissa bit    */
  221. }
  222.  
  223. long rnd_i ()
  224. /*
  225.  * returns a positive, uniformly distributed random number in [0,0x7fffffff]
  226.  */
  227. { 
  228.     register unsigned long i, j, *t = Fib;
  229.  
  230.     i = Fib_ind;
  231.     j = t[i];                    /* = X(n-55) */
  232.     j -= (i >= 24) ? t[i - 24] : t[i + 21]; /* = X(n-24) */
  233.     t[i] = j;
  234.     if (++i >= 55) i = 0;
  235.     Fib_ind = i;
  236.  
  237.     t = &auxtab[(j >> 24) & 0xff];
  238.     i = *t;
  239.     Xrnd_var = *t = Xrnd(Xrnd_var);
  240.     t = &prmtab[j & 0x3c];
  241.  
  242.     j =  *t++ & i;
  243.     j |= *t++ & ((i << 24) | ((i >>  8) & 0x00ffffff));
  244.     j |= *t++ & ((i << 16) | ((i >> 16) & 0x0000ffff));
  245.     j |= *t   & ((i <<  8) | ((i >> 24) & 0x000000ff));
  246.     
  247.     return j & MAX_INT;
  248. }
  249.  
  250.  
  251. long rnd_ri (rng)
  252.     long rng;
  253. /*
  254.  * randint: Return a random integer in a given Range [0..rng-1]
  255.  *          Note:  0 < rng
  256.  */
  257. {
  258.     register unsigned long  r, a;
  259.  
  260.     do {
  261.     r = rnd_i();
  262.     a = (r / rng) + 1;
  263.     a *= rng;
  264.     } while (a >= MAX_INT);
  265.     
  266.     a--;
  267.     return a - r;
  268. }
  269.