share/bobrand.c

   1 /* Author: Sam Trahan, October 2011
   2
   3    This is the implementation of a good but simple random number
   4    generator designed by Bob Jenkins, who placed it in the public
   5    domain.  His website described the algorithm and its public domain
   6    status on 2:23 AM EDT October 8, 2011 at this location:
   7
   8        http://burtleburtle.net/bob/rand/smallprng.html
   9
  10    And at that time, it said, "I wrote this PRNG. I place it in the
  11    public domain."  (PNRG is an acronym for "psuedo-random number
  12    generator" as defined elsewhere on his website.)
  13
  14    I modified his code to work as an array of random number generators
  15    and generate four output types (float, double, int32, int64).  This
  16    code is tested on the Intel, IBM and GNU C compilers, and will
  17    successfully produce identical floating-point numbers in [0,1) on
  18    all three compilers.  This code is not sensitive to optimization
  19    since all calculations are integer calculations, and hence are
  20    exact.
  21
  22    This algorithm, unlike the common Mersenne Twister, is not
  23    cryptographically secure, so don't use it to encrypt your banking
  24    information.  However, it does pass the entire suite of DIEHARD
  25    tests, so it is sufficiently random for meterological purposes.
  26    Its advantage over cryptographically secure algorithms is that it
  27    only needs 16 bytes to store its state, and is very fast, allowing
  28    us to have an independent random number generator for each
  29    gridpoint.  That avoids domain decomposition issues and allows us
  30    to generate random numbers in parallel across all processes,
  31    producing the same results regardless of which process or thread
  32    has which gridpoint.
  33
  34    Don't change any of the constants in this file without rerunning
  35    the full suite of randomness tests as described on Bob's website.
  36    Also, don't change the floating-point conversion unless you first
  37    test that it correctly produces 0, never produces 1.0, is uniformly
  38    distributed, and produces identical results on at least the Intel,
  39    GNU and IBM C compilers.
  40 */
  41 #include <stdint.h>
  42
  43 typedef uint32_t u4;
  44 typedef uint64_t u8;
  45
  46 #define rot(x,k) (((x)<<(k))|((x)>>(32-(k))))
  47
  48 void bobranval_impl( u4 *a, u4 *b, u4 *c, u4 *d, u4 *n ) {
  49   u4 e,i,nd=*n;
  50
  51   for(i=0;i<nd;i++) {
  52     e = a[i] - rot(b[i], 27);
  53     a[i] = b[i] ^ rot(c[i], 17);
  54     b[i] = c[i] + d[i];
  55     c[i] = d[i] + e;
  56     d[i] = e + a[i];
  57   }
  58 }
  59
  60 void bob_int_hash(u4 *in, u4 *out) {
  61   u4 a=0xf1ea5eed;
  62   u4 b,c,d,e,i;
  63   b=c=d=*in;
  64
  65   for(i=0;i<20;i++) {
  66     e = a - rot(b, 27);
  67     a = b ^ rot(c, 17);
  68     b = c + d;
  69     c = d + e;
  70     d = e + a;
  71   }
  72
  73   *out=d;
  74 }
  75
  76 void bobranval_r4_impl( u4 *a, u4 *b, u4 *c, u4 *d, float *result, u4 *n ) {
  77   /* 32-bit floating point implementation */
  78   u4 i,nd=*n;
  79
  80   bobranval_impl(a,b,c,d,n);
  81   for(i=0;i<nd;i++) {
  82     result[i]=(d[i]&0xfffff000)*2.328305e-10f;
  83   }
  84 }
  85
  86 void bobranval_i4_impl( u4 *a, u4 *b, u4 *c, u4 *d, u4 *result, u4 *n ) {
  87   /* 32-bit integer implementation */
  88   u4 i,nd=*n;
  89
  90   bobranval_impl(a,b,c,d,n);
  91   for(i=0;i<nd;i++)
  92     result[i]=d[i];
  93 }
  94
  95 void bobranval_i8_impl( u4 *a, u4 *b, u4 *c, u4 *d, u8 *result, u4 *n ) {
  96   /* 64-bit integer implementation */
  97   u4 i,nd=*n;
  98
  99   bobranval_impl(a,b,c,d,n);
 100   for(i=0;i<nd;i++)
 101     result[i]=d[i];
 102
 103   bobranval_impl(a,b,c,d,n);
 104   for(i=0;i<nd;i++)
 105     result[i]=(result[i]<<32) | d[i];
 106 }
 107
 108 void bobranval_r8_impl( u4 *a, u4 *b, u4 *c, u4 *d, u8 *result, u4 *n ) {
 109   /* 64-bit floating-point implementation */
 110   u4 i,nd=*n;
 111
 112   bobranval_impl(a,b,c,d,n);
 113   for(i=0;i<nd;i++)
 114     result[i]=d[i];
 115
 116   bobranval_impl(a,b,c,d,n);
 117   for(i=0;i<nd;i++) {
 118     ((double*)result)[i] =
 119        (((result[i]<<32) | d[i]) & 0xfffffffffffffc00ll)
 120              * 5.4210108624275218691107101938441e-20;
 121   }
 122 }
 123
 124 void bobraninit(u4 *a, u4 *b, u4 *c, u4 *d, u4 *seeds, u4 *seed2, u4 *n) {
 125   u4 i,nd=*n,one=1,iter;
 126   for(i=0;i<nd;i++) {
 127     a[i] = 0xf1ea5eed;
 128     b[i] = c[i] = d[i] = seeds[i]^*seed2;
 129     for (iter=0; iter<20; ++iter) {
 130       bobranval_impl(a+i,b+i,c+i,d+i,&one);
 131     }
 132   }
 133 }
 134
 135 /* Aliases for various fortran compilers */
 136
 137 void int_hash(u4*i,u4*o) { bob_int_hash(i,o); }
 138 void int_hash_(u4*i,u4*o) { bob_int_hash(i,o); }
 139 void int_hash__(u4*i,u4*o) { bob_int_hash(i,o); }
 140 void INT_HASH(u4*i,u4*o) { bob_int_hash(i,o); }
 141 void INT_HASH_(u4*i,u4*o) { bob_int_hash(i,o); }
 142 void INT_HASH__(u4*i,u4*o) { bob_int_hash(i,o); }
 143
 144 void bobraninit_(u4*a,u4*b,u4*c,u4*d,u4*s,u4*s2,u4*n) { bobraninit(a,b,c,d,s,s2,n); }
 145 void bobraninit__(u4*a,u4*b,u4*c,u4*d,u4*s,u4*s2,u4*n) { bobraninit(a,b,c,d,s,s2,n); }
 146 void BOBRANINIT_(u4*a,u4*b,u4*c,u4*d,u4*s,u4*s2,u4*n) { bobraninit(a,b,c,d,s,s2,n); }
 147 void BOBRANINIT__(u4*a,u4*b,u4*c,u4*d,u4*s,u4*s2,u4*n) { bobraninit(a,b,c,d,s,s2,n); }
 148
 149 void bobranval_r4(u4*a,u4*b,u4*c,u4*d,float*f,u4*n) { bobranval_r4_impl(a,b,c,d,f,n); }
 150 void bobranval_r4_(u4*a,u4*b,u4*c,u4*d,float*f,u4*n) { bobranval_r4_impl(a,b,c,d,f,n); }
 151 void bobranval_r4__(u4*a,u4*b,u4*c,u4*d,float*f,u4*n) { bobranval_r4_impl(a,b,c,d,f,n); }
 152 void BOBRANVAL_R4_(u4*a,u4*b,u4*c,u4*d,float*f,u4*n) { bobranval_r4_impl(a,b,c,d,f,n); }
 153 void BOBRANVAL_R4__(u4*a,u4*b,u4*c,u4*d,float*f,u4*n) { bobranval_r4_impl(a,b,c,d,f,n); }
 154
 155 void bobranval_i4(u4*a,u4*b,u4*c,u4*d,u4*i,u4*n) { bobranval_i4_impl(a,b,c,d,i,n); }
 156 void bobranval_i4_(u4*a,u4*b,u4*c,u4*d,u4*i,u4*n) { bobranval_i4_impl(a,b,c,d,i,n); }
 157 void bobranval_i4__(u4*a,u4*b,u4*c,u4*d,u4*i,u4*n) { bobranval_i4_impl(a,b,c,d,i,n); }
 158 void BOBRANVAL_I4_(u4*a,u4*b,u4*c,u4*d,u4*i,u4*n) { bobranval_i4_impl(a,b,c,d,i,n); }
 159 void BOBRANVAL_I4__(u4*a,u4*b,u4*c,u4*d,u4*i,u4*n) { bobranval_i4_impl(a,b,c,d,i,n); }
 160
 161 void bobranval_r8(u4*a,u4*b,u4*c,u4*d,u8*f,u4*n) { bobranval_r8_impl(a,b,c,d,f,n); }
 162 void bobranval_r8_(u4*a,u4*b,u4*c,u4*d,u8*f,u4*n) { bobranval_r8_impl(a,b,c,d,f,n); }
 163 void bobranval_r8__(u4*a,u4*b,u4*c,u4*d,u8*f,u4*n) { bobranval_r8_impl(a,b,c,d,f,n); }
 164 void BOBRANVAL_R8_(u4*a,u4*b,u4*c,u4*d,u8*f,u4*n) { bobranval_r8_impl(a,b,c,d,f,n); }
 165 void BOBRANVAL_R8__(u4*a,u4*b,u4*c,u4*d,u8*f,u4*n) { bobranval_r8_impl(a,b,c,d,f,n); }
 166
 167 void bobranval_i8(u4*a,u4*b,u4*c,u4*d,u8*i,u4*n) { bobranval_i8_impl(a,b,c,d,i,n); }
 168 void bobranval_i8_(u4*a,u4*b,u4*c,u4*d,u8*i,u4*n) { bobranval_i8_impl(a,b,c,d,i,n); }
 169 void bobranval_i8__(u4*a,u4*b,u4*c,u4*d,u8*i,u4*n) { bobranval_i8_impl(a,b,c,d,i,n); }
 170 void BOBRANVAL_I8_(u4*a,u4*b,u4*c,u4*d,u8*i,u4*n) { bobranval_i8_impl(a,b,c,d,i,n); }
 171 void BOBRANVAL_I8__(u4*a,u4*b,u4*c,u4*d,u8*i,u4*n) { bobranval_i8_impl(a,b,c,d,i,n); }