git.sesse.net Git - wloh/blob - 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 }