9 # include "ziggurat.hpp"
11 //****************************************************************************80
13 float r4_exp ( unsigned long int *jsr, int ke[256], float fe[256],
16 //****************************************************************************80
20 // R4_EXP returns an exponentially distributed single precision real value.
24 // The underlying algorithm is the ziggurat method.
26 // Before the first call to this function, the user must call R4_EXP_SETUP
27 // to determine the values of KE, FE and WE.
31 // This code is distributed under the GNU LGPL license.
43 // George Marsaglia, Wai Wan Tsang,
44 // The Ziggurat Method for Generating Random Variables,
45 // Journal of Statistical Software,
46 // Volume 5, Number 8, October 2000, seven pages.
50 // Input/output, unsigned long int *JSR, the seed.
52 // Input, int KE[256], data computed by R4_EXP_SETUP.
54 // Input, float FE[256], WE[256], data computed by R4_EXP_SETUP.
56 // Output, float R4_EXP, an exponentially distributed random value.
67 if ( abs ( jz ) < ke[iz] )
69 value = ( float ) ( abs ( jz ) ) * we[iz];
77 value = 7.69711 - log ( r4_uni ( jsr ) );
81 x = ( float ) ( abs ( jz ) ) * we[iz];
83 if ( fe[iz] + r4_uni ( jsr ) * ( fe[iz-1] - fe[iz] ) < exp ( - x ) )
92 if ( abs ( jz ) < ke[iz] )
94 value = ( float ) ( abs ( jz ) ) * we[iz];
101 //****************************************************************************80
103 void r4_exp_setup ( int ke[256], float fe[256], float we[256] )
105 //****************************************************************************80
109 // R4_EXP_SETUP sets data needed by R4_EXP.
113 // This code is distributed under the GNU LGPL license.
125 // George Marsaglia, Wai Wan Tsang,
126 // The Ziggurat Method for Generating Random Variables,
127 // Journal of Statistical Software,
128 // Volume 5, Number 8, October 2000, seven pages.
132 // Output, int KE[256], data needed by R4_EXP.
134 // Output, float FE[256], WE[256], data needed by R4_EXP.
137 double de = 7.697117470131487;
139 const double m2 = 2147483648.0;
141 double te = 7.697117470131487;
142 const double ve = 3.949659822581572E-03;
144 q = ve / exp ( - de );
146 ke[0] = ( int ) ( ( de / q ) * m2 );
149 we[0] = ( float ) ( q / m2 );
150 we[255] = ( float ) ( de / m2 );
153 fe[255] = ( float ) ( exp ( - de ) );
155 for ( i = 254; 1 <= i; i-- )
157 de = - log ( ve / de + exp ( - de ) );
158 ke[i+1] = ( int ) ( ( de / te ) * m2 );
160 fe[i] = ( float ) ( exp ( - de ) );
161 we[i] = ( float ) ( de / m2 );
165 //****************************************************************************80
167 float r4_nor ( unsigned long int *jsr, int kn[128], float fn[128],
170 //****************************************************************************80
174 // R4_NOR returns a normally distributed single precision real value.
178 // The value returned is generated from a distribution with mean 0 and
181 // The underlying algorithm is the ziggurat method.
183 // Before the first call to this function, the user must call R4_NOR_SETUP
184 // to determine the values of KN, FN and WN.
188 // This code is distributed under the GNU LGPL license.
200 // George Marsaglia, Wai Wan Tsang,
201 // The Ziggurat Method for Generating Random Variables,
202 // Journal of Statistical Software,
203 // Volume 5, Number 8, October 2000, seven pages.
207 // Input/output, unsigned long int *JSR, the seed.
209 // Input, int KN[128], data computed by R4_NOR_SETUP.
211 // Input, float FN[128], WN[128], data computed by R4_NOR_SETUP.
213 // Output, float R4_NOR, a normally distributed random value.
218 const float r = 3.442620;
226 if ( abs ( hz ) < kn[iz] )
228 value = ( float ) ( hz ) * wn[iz];
238 x = - 0.2904764 * log ( r4_uni ( jsr ) );
239 y = - log ( r4_uni ( jsr ) );
240 if ( x * x <= y + y );
257 x = ( float ) ( hz ) * wn[iz];
259 if ( fn[iz] + r4_uni ( jsr ) * ( fn[iz-1] - fn[iz] ) < exp ( - 0.5 * x * x ) )
268 if ( abs ( hz ) < kn[iz] )
270 value = ( float ) ( hz ) * wn[iz];
278 //****************************************************************************80
280 void r4_nor_setup ( int kn[128], float fn[128], float wn[128] )
282 //****************************************************************************80
286 // R4_NOR_SETUP sets data needed by R4_NOR.
290 // This code is distributed under the GNU LGPL license.
302 // George Marsaglia, Wai Wan Tsang,
303 // The Ziggurat Method for Generating Random Variables,
304 // Journal of Statistical Software,
305 // Volume 5, Number 8, October 2000, seven pages.
309 // Output, int KN[128], data needed by R4_NOR.
311 // Output, float FN[128], WN[128], data needed by R4_NOR.
314 double dn = 3.442619855899;
316 const double m1 = 2147483648.0;
318 double tn = 3.442619855899;
319 const double vn = 9.91256303526217E-03;
321 q = vn / exp ( - 0.5 * dn * dn );
323 kn[0] = ( int ) ( ( dn / q ) * m1 );
326 wn[0] = ( float ) ( q / m1 );
327 wn[127] = ( float ) ( dn / m1 );
330 fn[127] = ( float ) ( exp ( - 0.5 * dn * dn ) );
332 for ( i = 126; 1 <= i; i-- )
334 dn = sqrt ( - 2.0 * log ( vn / dn + exp ( - 0.5 * dn * dn ) ) );
335 kn[i+1] = ( int ) ( ( dn / tn ) * m1 );
337 fn[i] = ( float ) ( exp ( - 0.5 * dn * dn ) );
338 wn[i] = ( float ) ( dn / m1 );
342 //****************************************************************************80
344 float r4_uni ( unsigned long int *jsr )
346 //****************************************************************************80
350 // R4_UNI returns a uniformly distributed real value.
354 // This code is distributed under the GNU LGPL license.
366 // George Marsaglia, Wai Wan Tsang,
367 // The Ziggurat Method for Generating Random Variables,
368 // Journal of Statistical Software,
369 // Volume 5, Number 8, October 2000, seven pages.
373 // Input/output, unsigned long int *JSR, the seed.
375 // Output, float R4_UNI, a uniformly distributed random value in
379 unsigned long int jsr_input;
384 *jsr = ( *jsr ^ ( *jsr << 13 ) );
385 *jsr = ( *jsr ^ ( *jsr >> 17 ) );
386 *jsr = ( *jsr ^ ( *jsr << 5 ) );
388 value = fmod ( 0.5 + ( float ) ( jsr_input + *jsr ) / 65536.0 / 65536.0, 1.0 );
392 //****************************************************************************80
394 unsigned long int shr3 ( unsigned long int *jsr )
396 //****************************************************************************80
400 // SHR3 evaluates the SHR3 generator for integers.
404 // This code is distributed under the GNU LGPL license.
416 // George Marsaglia, Wai Wan Tsang,
417 // The Ziggurat Method for Generating Random Variables,
418 // Journal of Statistical Software,
419 // Volume 5, Number 8, October 2000, seven pages.
423 // Input/output, unsigned long int *JSR, the seed, which is updated
426 // Output, unsigned long int SHR3, the new value.
429 unsigned long int value;
433 *jsr = ( *jsr ^ ( *jsr << 13 ) );
434 *jsr = ( *jsr ^ ( *jsr >> 17 ) );
435 *jsr = ( *jsr ^ ( *jsr << 5 ) );
437 value = value + *jsr;
441 //****************************************************************************80
445 //****************************************************************************80
449 // TIMESTAMP prints the current YMDHMS date as a time stamp.
453 // 31 May 2001 09:45:54 AM
457 // This code is distributed under the GNU LGPL license.
472 # define TIME_SIZE 40
474 static char time_buffer[TIME_SIZE];
480 tm = localtime ( &now );
482 len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
484 cout << time_buffer << "\n";