27 #define UINT64_LITERAL(X) X ## ui64 29 #define UINT64_LITERAL(X) X ## ULL 32 #define RAND48_SEED UINT64_LITERAL(0x00001234abcd330e) 33 #define RAND48_MULT UINT64_LITERAL(0x00000005deece66d) 34 #define RAND48_ADD UINT64_LITERAL(0x000000000000000b) 35 #define RAND48_MASK UINT64_LITERAL(0x0000ffffffffffff) 41 double second_gaussian;
42 uint64_t second_gaussian_waiting;
61 void init(
unsigned long seed) {
63 second_gaussian_waiting = 0;
65 rand48_seed = rand48_seed << 16;
73 rand48_seed = ( rand48_seed * rand48_mult + rand48_add ) &
RAND48_MASK;
77 void split(
int iStream,
int numStreams) {
85 for ( i = 0; i < iStream; ++i )
skip();
88 uint64_t save_seed = rand48_seed;
91 rand48_seed = rand48_add;
92 for ( i = 1; i < numStreams; ++i )
skip();
93 uint64_t new_add = rand48_seed;
96 rand48_seed = rand48_mult;
98 for ( i = 1; i < numStreams; ++i )
skip();
99 rand48_mult = rand48_seed;
101 rand48_add = new_add;
102 rand48_seed = save_seed;
105 second_gaussian_waiting = 0;
112 return ( (
double) rand48_seed * exp48 );
118 static int count = 0;
122 if (second_gaussian_waiting) {
123 second_gaussian_waiting = 0;
126 fprintf(stderr,
"RANDOM %6d %12.9f\n", count, second_gaussian);
128 return second_gaussian;
131 while (r >=1. || r < 1.523e-8) {
136 fac = sqrt(-2.0 * log(r)/r);
139 second_gaussian_waiting = 1;
140 second_gaussian = v1 * fac;
143 fprintf(stderr,
"RANDOM %6d %12.9f\n", count, v2 * fac);
192 if (num_gaussians <= 2) {
194 for(
int i=0; i<num_gaussians; ++i) {
199 }
else if (2 < num_gaussians && num_gaussians <= 200) {
200 const BigReal d = 0.5*num_gaussians - 1./3;
201 const BigReal c = 1. / (3*sqrt(d));
209 v = (1 + c*z); v = v*v*v;
211 while(log(u) >= (0.5*z*z + d*(1 - v + log(v))));
214 return num_gaussians + sqrt(2*num_gaussians)*
gaussian();
227 return ( ( rand48_seed >> 17 ) &
UINT64_LITERAL(0x000000007fffffff) );
231 template <
class Elem>
void reorder(Elem *a,
int n) {
232 for (
int i = 0; i < (n-1); ++i ) {
233 int ie = i + (
integer() % (n-i) );
234 if ( ie == i )
continue;
235 const Elem e = a[ie];
247 for (
int i=0; i < n; i++) {
258 for (
int i=0; i < n; i++) {
Vector gaussian_vector(void)
BigReal sum_of_squared_gaussians(int64_t num_gaussians)
Return the sum of i.i.d. squared Gaussians.
void uniform_array_f(float *a, int n)
void split(int iStream, int numStreams)
void reorder(Elem *a, int n)
void init(unsigned long seed)
Random(unsigned long seed)
void gaussian_array_f(float *a, int n)
#define UINT64_LITERAL(X)