]> git.sesse.net Git - wloh/blob - ziggurat.cpp
Show rating standard deviation.
[wloh] / ziggurat.cpp
1 # include <cstdlib>
2 # include <iostream>
3 # include <iomanip>
4 # include <cmath>
5 # include <ctime>
6
7 using namespace std;
8
9 # include "ziggurat.hpp"
10
11 //****************************************************************************80
12
13 float r4_exp ( unsigned long int *jsr, int ke[256], float fe[256], 
14   float we[256] )
15
16 //****************************************************************************80
17 //
18 //  Purpose:
19 //
20 //    R4_EXP returns an exponentially distributed single precision real value.
21 //
22 //  Discussion:
23 //
24 //    The underlying algorithm is the ziggurat method.
25 //
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.
28 //
29 //  Licensing:
30 //
31 //    This code is distributed under the GNU LGPL license. 
32 //
33 //  Modified:
34 //
35 //    08 December 20080
36 //
37 //  Author:
38 //
39 //    John Burkardt
40 //
41 //  Reference:
42 //
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.
47 //
48 //  Parameters:
49 //
50 //    Input/output, unsigned long int *JSR, the seed.
51 //
52 //    Input, int KE[256], data computed by R4_EXP_SETUP.
53 //
54 //    Input, float FE[256], WE[256], data computed by R4_EXP_SETUP.
55 //
56 //    Output, float R4_EXP, an exponentially distributed random value.
57 //
58 {
59   int iz;
60   int jz;
61   float value;
62   float x;
63
64   jz = shr3 ( jsr );
65   iz = ( jz & 255 );
66
67   if ( abs ( jz  ) < ke[iz] )
68   {
69     value = ( float ) ( abs ( jz ) ) * we[iz];
70   }
71   else
72   {
73     for ( ; ; )
74     {
75       if ( iz == 0 )
76       {
77         value = 7.69711 - log ( r4_uni ( jsr ) );
78         break;
79       }
80
81       x = ( float ) ( abs ( jz ) ) * we[iz];
82
83       if ( fe[iz] + r4_uni ( jsr ) * ( fe[iz-1] - fe[iz] ) < exp ( - x ) )
84       {
85         value = x;
86         break;
87       }
88
89       jz = shr3 ( jsr );
90       iz = ( jz & 255 );
91
92       if ( abs ( jz ) < ke[iz] )
93       {
94         value = ( float ) ( abs ( jz ) ) * we[iz];
95         break;
96       }
97     }
98   }
99   return value;
100 }
101 //****************************************************************************80
102
103 void r4_exp_setup ( int ke[256], float fe[256], float we[256] )
104
105 //****************************************************************************80
106 //
107 //  Purpose:
108 //
109 //    R4_EXP_SETUP sets data needed by R4_EXP.
110 //
111 //  Licensing:
112 //
113 //    This code is distributed under the GNU LGPL license. 
114 //
115 //  Modified:
116 //
117 //    08 December 2008
118 //
119 //  Author:
120 //
121 //    John Burkardt
122 //
123 //  Reference:
124 //
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.
129 //
130 //  Parameters:
131 //
132 //    Output, int KE[256], data needed by R4_EXP.
133 //
134 //    Output, float FE[256], WE[256], data needed by R4_EXP.
135 //
136 {
137   double de = 7.697117470131487;
138   int i;
139   const double m2 = 2147483648.0;
140   double q;
141   double te = 7.697117470131487;
142   const double ve = 3.949659822581572E-03;
143
144   q = ve / exp ( - de );
145
146   ke[0] = ( int ) ( ( de / q ) * m2 );
147   ke[1] = 0;
148
149   we[0] = ( float ) ( q / m2 );
150   we[255] = ( float ) ( de / m2 );
151
152   fe[0] = 1.0;
153   fe[255] = ( float ) ( exp ( - de ) );
154
155   for ( i = 254; 1 <= i; i-- )
156   {
157     de = - log ( ve / de + exp ( - de ) );
158     ke[i+1] = ( int ) ( ( de / te ) * m2 );
159     te = de;
160     fe[i] = ( float ) ( exp ( - de ) );
161     we[i] = ( float ) ( de / m2 );
162   }
163   return;
164 }
165 //****************************************************************************80
166
167 float r4_nor ( unsigned long int *jsr, int kn[128], float fn[128], 
168   float wn[128] )
169
170 //****************************************************************************80
171 //
172 //  Purpose:
173 //
174 //    R4_NOR returns a normally distributed single precision real value.
175 //
176 //  Discussion:
177 //
178 //    The value returned is generated from a distribution with mean 0 and 
179 //    variance 1.
180 //
181 //    The underlying algorithm is the ziggurat method.
182 //
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.
185 //
186 //  Licensing:
187 //
188 //    This code is distributed under the GNU LGPL license. 
189 //
190 //  Modified:
191 //
192 //    08 December 2008
193 //
194 //  Author:
195 //
196 //    John Burkardt
197 //
198 //  Reference:
199 //
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.
204 //
205 //  Parameters:
206 //
207 //    Input/output, unsigned long int *JSR, the seed.
208 //
209 //    Input, int KN[128], data computed by R4_NOR_SETUP.
210 //
211 //    Input, float FN[128], WN[128], data computed by R4_NOR_SETUP.
212 //
213 //    Output, float R4_NOR, a normally distributed random value.
214 //
215 {
216   int hz;
217   int iz;
218   const float r = 3.442620;
219   float value;
220   float x;
221   float y;
222
223   hz = shr3 ( jsr );
224   iz = ( hz & 127 );
225
226   if ( abs ( hz ) < kn[iz] )
227   {
228     value = ( float ) ( hz ) * wn[iz];
229   }
230   else
231   {
232     for ( ; ; )
233     {
234       if ( iz == 0 )
235       {
236         for ( ; ; )
237         {
238           x = - 0.2904764 * log ( r4_uni ( jsr ) );
239           y = - log ( r4_uni ( jsr ) );
240           if ( x * x <= y + y );
241           {
242             break;
243           }
244         }
245
246         if ( hz <= 0 )
247         {
248           value = - r - x;
249         }
250         else
251         {
252           value = + r + x;
253         }
254         break;
255       }
256
257       x = ( float ) ( hz ) * wn[iz];
258
259       if ( fn[iz] + r4_uni ( jsr ) * ( fn[iz-1] - fn[iz] ) < exp ( - 0.5 * x * x ) )
260       {
261         value = x;
262         break;
263       }
264
265       hz = shr3 ( jsr );
266       iz = ( hz & 127 );
267
268       if ( abs ( hz ) < kn[iz] )
269       {
270         value = ( float ) ( hz ) * wn[iz];
271         break;
272       }
273     }
274   }
275
276   return value;
277 }
278 //****************************************************************************80
279
280 void r4_nor_setup ( int kn[128], float fn[128], float wn[128] )
281
282 //****************************************************************************80
283 //
284 //  Purpose:
285 //
286 //    R4_NOR_SETUP sets data needed by R4_NOR.
287 //
288 //  Licensing:
289 //
290 //    This code is distributed under the GNU LGPL license. 
291 //
292 //  Modified:
293 //
294 //    04 May 2008
295 //
296 //  Author:
297 //
298 //    John Burkardt
299 //
300 //  Reference:
301 //
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.
306 //
307 //  Parameters:
308 //
309 //    Output, int KN[128], data needed by R4_NOR.
310 //
311 //    Output, float FN[128], WN[128], data needed by R4_NOR.
312 //
313 {
314   double dn = 3.442619855899;
315   int i;
316   const double m1 = 2147483648.0;
317   double q;
318   double tn = 3.442619855899;
319   const double vn = 9.91256303526217E-03;
320
321   q = vn / exp ( - 0.5 * dn * dn );
322
323   kn[0] = ( int ) ( ( dn / q ) * m1 );
324   kn[1] = 0;
325
326   wn[0] = ( float ) ( q / m1 );
327   wn[127] = ( float ) ( dn / m1 );
328
329   fn[0] = 1.0;
330   fn[127] = ( float ) ( exp ( - 0.5 * dn * dn ) );
331
332   for ( i = 126; 1 <= i; i-- )
333   {
334     dn = sqrt ( - 2.0 * log ( vn / dn + exp ( - 0.5 * dn * dn ) ) );
335     kn[i+1] = ( int ) ( ( dn / tn ) * m1 );
336     tn = dn;
337     fn[i] = ( float ) ( exp ( - 0.5 * dn * dn ) );
338     wn[i] = ( float ) ( dn / m1 );
339   }
340   return;
341 }
342 //****************************************************************************80
343
344 float r4_uni ( unsigned long int *jsr )
345
346 //****************************************************************************80
347 //
348 //  Purpose:
349 //
350 //    R4_UNI returns a uniformly distributed real value.
351 //
352 //  Licensing:
353 //
354 //    This code is distributed under the GNU LGPL license. 
355 //
356 //  Modified:
357 //
358 //    20 May 2008
359 //
360 //  Author:
361 //
362 //    John Burkardt
363 //
364 //  Reference:
365 //
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.
370 //
371 //  Parameters:
372 //
373 //    Input/output, unsigned long int *JSR, the seed.
374 //
375 //    Output, float R4_UNI, a uniformly distributed random value in
376 //    the range [0,1].
377 //
378 {
379   unsigned long int jsr_input;
380   float value;
381
382   jsr_input = *jsr;
383
384   *jsr = ( *jsr ^ ( *jsr <<   13 ) );
385   *jsr = ( *jsr ^ ( *jsr >>   17 ) );
386   *jsr = ( *jsr ^ ( *jsr <<    5 ) );
387
388   value = fmod ( 0.5 + ( float ) ( jsr_input + *jsr ) / 65536.0 / 65536.0, 1.0 );
389
390   return value;
391 }
392 //****************************************************************************80
393
394 unsigned long int shr3 ( unsigned long int *jsr )
395
396 //****************************************************************************80
397 //
398 //  Purpose:
399 //
400 //    SHR3 evaluates the SHR3 generator for integers.
401 //
402 //  Licensing:
403 //
404 //    This code is distributed under the GNU LGPL license. 
405 //
406 //  Modified:
407 //
408 //    08 December 2008
409 //
410 //  Author:
411 //
412 //    John Burkardt
413 //
414 //  Reference:
415 //
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.
420 //
421 //  Parameters:
422 //
423 //    Input/output, unsigned long int *JSR, the seed, which is updated 
424 //    on each call.
425 //
426 //    Output, unsigned long int SHR3, the new value.
427 //
428 {
429   unsigned long int value;
430
431   value = *jsr;
432
433   *jsr = ( *jsr ^ ( *jsr <<   13 ) );
434   *jsr = ( *jsr ^ ( *jsr >>   17 ) );
435   *jsr = ( *jsr ^ ( *jsr <<    5 ) );
436
437   value = value + *jsr;
438
439   return value;
440 }
441 //****************************************************************************80
442
443 void timestamp ( )
444
445 //****************************************************************************80
446 //
447 //  Purpose:
448 //
449 //    TIMESTAMP prints the current YMDHMS date as a time stamp.
450 //
451 //  Example:
452 //
453 //    31 May 2001 09:45:54 AM
454 //
455 //  Licensing:
456 //
457 //    This code is distributed under the GNU LGPL license. 
458 //
459 //  Modified:
460 //
461 //    24 September 2003
462 //
463 //  Author:
464 //
465 //    John Burkardt
466 //
467 //  Parameters:
468 //
469 //    None
470 //
471 {
472 # define TIME_SIZE 40
473
474   static char time_buffer[TIME_SIZE];
475   const struct tm *tm;
476   size_t len;
477   time_t now;
478
479   now = time ( NULL );
480   tm = localtime ( &now );
481
482   len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
483
484   cout << time_buffer << "\n";
485
486   return;
487 # undef TIME_SIZE
488 }