limit fstBC to 30bp in Python3 ver.
[GalaxyCodeBases.git] / BGI / soap_src / soap_builder / r250.c
blob55311554feff2385d847663841bf1afaddf9e3bd
1 /* r250.c the r250 uniform random number algorithm
3 Kirkpatrick, S., and E. Stoll, 1981; "A Very Fast
4 Shift-Register Sequence Random Number Generator",
5 Journal of Computational Physics, V.40
7 also:
9 see W.L. Maier, DDJ May 1991
15 #include <limits.h>
16 #include "r250.h"
18 // static functions
19 static unsigned int randlcg();
22 #define MSB 0x80000000L
23 #define ALL_BITS 0xffffffffL
24 #define HALF_RANGE 0x40000000L
25 #define STEP 7
26 #define BITS 32
28 static unsigned int r250_buffer[ 250 ];
29 static int r250_index;
31 static unsigned int randlcg(int sd) /* returns a random unsigned integer */
33 static long quotient1 = LONG_MAX / 16807L;
34 static long remainder1 = LONG_MAX % 16807L;
36 if ( sd <= quotient1 )
37 sd = (sd * 16807L) % LONG_MAX;
38 else
40 int high_part = sd / quotient1;
41 int low_part = sd % quotient1;
43 int test = 16807L * low_part - remainder1 * high_part;
45 if ( test > 0 )
46 sd = test;
47 else
48 sd = test + LONG_MAX;
52 return sd;
55 void r250_init(int sd)
57 int j, k;
58 unsigned int mask, msb;
60 r250_index = 0;
61 for (j = 0; j < 250; j++) /* fill r250 buffer with BITS-1 bit values */
62 sd = r250_buffer[j] = randlcg(sd);
65 for (j = 0; j < 250; j++) /* set some MSBs to 1 */
66 if ( (sd=randlcg(sd)) > HALF_RANGE )
67 r250_buffer[j] |= MSB;
70 msb = MSB; /* turn on diagonal bit */
71 mask = ALL_BITS; /* turn off the leftmost bits */
73 for (j = 0; j < BITS; j++)
75 k = STEP * j + 3; /* select a word to operate on */
76 r250_buffer[k] &= mask; /* turn off bits left of the diagonal */
77 r250_buffer[k] |= msb; /* turn on the diagonal bit */
78 mask >>= 1;
79 msb >>= 1;
84 unsigned int r250() /* returns a random unsigned integer */
86 register int j;
87 register unsigned int new_rand;
89 if ( r250_index >= 147 )
90 j = r250_index - 147; /* wrap pointer around */
91 else
92 j = r250_index + 103;
94 new_rand = r250_buffer[ r250_index ] ^ r250_buffer[ j ];
95 r250_buffer[ r250_index ] = new_rand;
97 if ( r250_index >= 249 ) /* increment pointer for next time */
98 r250_index = 0;
99 else
100 r250_index++;
102 return new_rand;
107 double dr250() /* returns a random double in range 0..1 */
109 register int j;
110 register unsigned int new_rand;
112 if ( r250_index >= 147 )
113 j = r250_index - 147; /* wrap pointer around */
114 else
115 j = r250_index + 103;
117 new_rand = r250_buffer[ r250_index ] ^ r250_buffer[ j ];
118 r250_buffer[ r250_index ] = new_rand;
120 if ( r250_index >= 249 ) /* increment pointer for next time */
121 r250_index = 0;
122 else
123 r250_index++;
125 return (double)new_rand / ALL_BITS;