htswanalysis/MACS/lib/gsl/gsl-1.11/rng/test.c

   1 /* rng/test.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <stdlib.h>
  22 #include <stdio.h>
  23 #include <math.h>
  24 #include <gsl/gsl_rng.h>
  25 #include <gsl/gsl_test.h>
  26 #include <gsl/gsl_ieee_utils.h>
  27
  28 void rng_test (const gsl_rng_type * T, unsigned long int seed, unsigned int n,
  29                unsigned long int result);
  30 void rng_float_test (const gsl_rng_type * T);
  31 void generic_rng_test (const gsl_rng_type * T);
  32 void rng_state_test (const gsl_rng_type * T);
  33 void rng_parallel_state_test (const gsl_rng_type * T);
  34 void rng_read_write_test (const gsl_rng_type * T);
  35 int rng_max_test (gsl_rng * r, unsigned long int *kmax, unsigned long int ran_max) ;
  36 int rng_min_test (gsl_rng * r, unsigned long int *kmin, unsigned long int ran_min, unsigned long int ran_max) ;
  37 int rng_sum_test (gsl_rng * r, double *sigma);
  38 int rng_bin_test (gsl_rng * r, double *sigma);
  39
  40 #define N  10000
  41 #define N2 200000
  42
  43 int
  44 main (void)
  45 {
  46   const gsl_rng_type ** rngs = gsl_rng_types_setup();  /* get all rng types */
  47
  48   const gsl_rng_type ** r ;
  49
  50   gsl_ieee_env_setup ();
  51
  52   gsl_rng_env_setup ();
  53
  54   /* specific tests of known results for 10000 iterations with seed = 1 */
  55
  56   rng_test (gsl_rng_rand, 1, 10000, 1910041713);
  57   rng_test (gsl_rng_randu, 1, 10000, 1623524161);
  58   rng_test (gsl_rng_cmrg, 1, 10000, 719452880);
  59   rng_test (gsl_rng_minstd, 1, 10000, 1043618065);
  60   rng_test (gsl_rng_mrg, 1, 10000, 2064828650);
  61   rng_test (gsl_rng_taus, 1, 10000, 2733957125UL);
  62   rng_test (gsl_rng_taus2, 1, 10000, 2733957125UL);
  63   rng_test (gsl_rng_taus113, 1, 1000, 1925420673UL);
  64   rng_test (gsl_rng_transputer, 1, 10000, 1244127297UL);
  65   rng_test (gsl_rng_vax, 1, 10000, 3051034865UL);
  66
  67   /* Borosh13 test value from PARI: (1812433253^10000)%(2^32) */
  68   rng_test (gsl_rng_borosh13, 1, 10000, 2513433025UL);
  69
  70   /* Fishman18 test value from PARI: (62089911^10000)%(2^31-1) */
  71   rng_test (gsl_rng_fishman18, 1, 10000, 330402013UL);
  72
  73   /* Fishman2x test value from PARI:
  74      ((48271^10000)%(2^31-1) - (40692^10000)%(2^31-249))%(2^31-1) */
  75   rng_test (gsl_rng_fishman2x, 1, 10000, 540133597UL);
  76
  77   /* Knuthran2 test value from PARI:
  78      { xn1=1; xn2=1; for (n=1,10000,
  79             xn = (271828183*xn1 - 314159269*xn2)%(2^31-1);
  80             xn2=xn1; xn1=xn; print(xn); ) } */
  81   rng_test (gsl_rng_knuthran2, 1, 10000, 1084477620UL);
  82
  83   /* Knuthran test value taken from p188 in Knuth Vol 2. 3rd Ed */
  84   rng_test (gsl_rng_knuthran, 310952, 1009 * 2009 + 1, 461390032);
  85
  86   /* Knuthran improved test value from Knuth's source */
  87   rng_test (gsl_rng_knuthran2002, 310952, 1, 708622036);
  88   rng_test (gsl_rng_knuthran2002, 310952, 2, 1005450560);
  89   rng_test (gsl_rng_knuthran2002, 310952, 100 * 2009 + 1, 995235265);
  90   rng_test (gsl_rng_knuthran2002, 310952, 1009 * 2009 + 1, 704987132);
  91
  92   /* Lecuyer21 test value from PARI: (40692^10000)%(2^31-249) */
  93   rng_test (gsl_rng_lecuyer21, 1, 10000, 2006618587UL);
  94
  95   /* Waterman14 test value from PARI: (1566083941^10000)%(2^32) */
  96   rng_test (gsl_rng_waterman14, 1, 10000, 3776680385UL);
  97
  98   /* specific tests of known results for 10000 iterations with seed = 6 */
  99
 100   /* Coveyou test value from PARI:
 101      x=6; for(n=1,10000,x=(x*(x+1))%(2^32);print(x);) */
 102
 103   rng_test (gsl_rng_coveyou, 6, 10000, 1416754246UL);
 104
 105   /* Fishman20 test value from PARI: (6*48271^10000)%(2^31-1) */
 106   rng_test (gsl_rng_fishman20, 6, 10000, 248127575UL);
 107
 108   /* FIXME: the ranlux tests below were made by running the fortran code and
 109      getting the expected value from that. An analytic calculation
 110      would be preferable. */
 111
 112   rng_test (gsl_rng_ranlux, 314159265, 10000, 12077992);
 113   rng_test (gsl_rng_ranlux389, 314159265, 10000, 165942);
 114
 115   rng_test (gsl_rng_ranlxs0, 1, 10000, 11904320);
 116   /* 0.709552764892578125 * ldexp(1.0,24) */
 117
 118   rng_test (gsl_rng_ranlxs1, 1, 10000, 8734328);
 119   /* 0.520606517791748047 * ldexp(1.0,24) */
 120
 121   rng_test (gsl_rng_ranlxs2, 1, 10000, 6843140);
 122   /* 0.407882928848266602 * ldexp(1.0,24) */
 123
 124   rng_test (gsl_rng_ranlxd1, 1, 10000, 1998227290UL);
 125   /* 0.465248546261094020 * ldexp(1.0,32) */
 126
 127   rng_test (gsl_rng_ranlxd2, 1, 10000, 3949287736UL);
 128   /* 0.919515205581550532 * ldexp(1.0,32) */
 129
 130   /* FIXME: the tests below were made by running the original code in
 131      the ../random directory and getting the expected value from
 132      that. An analytic calculation would be preferable. */
 133
 134   rng_test (gsl_rng_slatec, 1, 10000, 45776);
 135   rng_test (gsl_rng_uni, 1, 10000, 9214);
 136   rng_test (gsl_rng_uni32, 1, 10000, 1155229825);
 137   rng_test (gsl_rng_zuf, 1, 10000, 3970);
 138
 139   /* The tests below were made by running the original code and
 140      getting the expected value from that. An analytic calculation
 141      would be preferable. */
 142
 143   rng_test (gsl_rng_r250, 1, 10000, 1100653588);
 144   rng_test (gsl_rng_mt19937, 4357, 1000, 1186927261);
 145   rng_test (gsl_rng_mt19937_1999, 4357, 1000, 1030650439);
 146   rng_test (gsl_rng_mt19937_1998, 4357, 1000, 1309179303);
 147   rng_test (gsl_rng_tt800, 0, 10000, 2856609219UL);
 148
 149   rng_test (gsl_rng_ran0, 0, 10000, 1115320064);
 150   rng_test (gsl_rng_ran1, 0, 10000, 1491066076);
 151   rng_test (gsl_rng_ran2, 0, 10000, 1701364455);
 152   rng_test (gsl_rng_ran3, 0, 10000, 186340785);
 153
 154   rng_test (gsl_rng_ranmar, 1, 10000, 14428370);
 155
 156   rng_test (gsl_rng_rand48, 0, 10000, 0xDE095043UL);
 157   rng_test (gsl_rng_rand48, 1, 10000, 0xEDA54977UL);
 158
 159   rng_test (gsl_rng_random_glibc2, 0, 10000, 1908609430);
 160   rng_test (gsl_rng_random8_glibc2, 0, 10000, 1910041713);
 161   rng_test (gsl_rng_random32_glibc2, 0, 10000, 1587395585);
 162   rng_test (gsl_rng_random64_glibc2, 0, 10000, 52848624);
 163   rng_test (gsl_rng_random128_glibc2, 0, 10000, 1908609430);
 164   rng_test (gsl_rng_random256_glibc2, 0, 10000, 179943260);
 165
 166   rng_test (gsl_rng_random_bsd, 0, 10000, 1457025928);
 167   rng_test (gsl_rng_random8_bsd, 0, 10000, 1910041713);
 168   rng_test (gsl_rng_random32_bsd, 0, 10000, 1663114331);
 169   rng_test (gsl_rng_random64_bsd, 0, 10000, 864469165);
 170   rng_test (gsl_rng_random128_bsd, 0, 10000, 1457025928);
 171   rng_test (gsl_rng_random256_bsd, 0, 10000, 1216357476);
 172
 173   rng_test (gsl_rng_random_libc5, 0, 10000, 428084942);
 174   rng_test (gsl_rng_random8_libc5, 0, 10000, 1910041713);
 175   rng_test (gsl_rng_random32_libc5, 0, 10000, 1967452027);
 176   rng_test (gsl_rng_random64_libc5, 0, 10000, 2106639801);
 177   rng_test (gsl_rng_random128_libc5, 0, 10000, 428084942);
 178   rng_test (gsl_rng_random256_libc5, 0, 10000, 116367984);
 179
 180   rng_test (gsl_rng_ranf, 0, 10000, 2152890433UL);
 181   rng_test (gsl_rng_ranf, 2, 10000, 339327233);
 182
 183   /* Test constant relationship between int and double functions */
 184
 185   for (r = rngs ; *r != 0; r++)
 186     rng_float_test (*r);
 187
 188   /* Test save/restore functions */
 189
 190   for (r = rngs ; *r != 0; r++)
 191     rng_state_test (*r);
 192
 193   for (r = rngs ; *r != 0; r++)
 194     rng_parallel_state_test (*r);
 195
 196   for (r = rngs ; *r != 0; r++)
 197     rng_read_write_test (*r);
 198
 199   /* generic statistical tests (these are just to make sure that we
 200      don't get any crazy results back from the generator, i.e. they
 201      aren't a test of the algorithm, just the implementation) */
 202
 203   for (r = rngs ; *r != 0; r++)
 204     generic_rng_test (*r);
 205
 206   exit (gsl_test_summary ());
 207 }
 208
 209
 210 void
 211 rng_test (const gsl_rng_type * T, unsigned long int seed, unsigned int n,
 212           unsigned long int result)
 213 {
 214   gsl_rng *r = gsl_rng_alloc (T);
 215   unsigned int i;
 216   unsigned long int k = 0;
 217   int status;
 218
 219   if (seed != 0)
 220     {
 221       gsl_rng_set (r, seed);
 222     }
 223
 224   for (i = 0; i < n; i++)
 225     {
 226       k = gsl_rng_get (r);
 227     }
 228
 229   status = (k != result);
 230   gsl_test (status, "%s, %u steps (%u observed vs %u expected)",
 231             gsl_rng_name (r), n, k, result);
 232
 233   gsl_rng_free (r);
 234 }
 235
 236 void
 237 rng_float_test (const gsl_rng_type * T)
 238 {
 239   gsl_rng *ri = gsl_rng_alloc (T);
 240   gsl_rng *rf = gsl_rng_alloc (T);
 241
 242   double u, c ;
 243   unsigned int i;
 244   unsigned long int k = 0;
 245   int status = 0 ;
 246
 247   do
 248     {
 249       k = gsl_rng_get (ri);
 250       u = gsl_rng_get (rf);
 251     }
 252   while (k == 0) ;
 253
 254   c = k / u ;
 255   for (i = 0; i < N2; i++)
 256     {
 257       k = gsl_rng_get (ri);
 258       u = gsl_rng_get (rf);
 259       if (c*k != u)
 260         {
 261           status = 1 ;
 262           break ;
 263         }
 264     }
 265
 266   gsl_test (status, "%s, ratio of int to double (%g observed vs %g expected)",
 267             gsl_rng_name (ri), c, k/u);
 268
 269   gsl_rng_free (ri);
 270   gsl_rng_free (rf);
 271 }
 272
 273
 274 void
 275 rng_state_test (const gsl_rng_type * T)
 276 {
 277   unsigned long int test_a[N], test_b[N];
 278
 279   int i;
 280
 281   gsl_rng *r = gsl_rng_alloc (T);
 282   gsl_rng *r_save = gsl_rng_alloc (T);
 283
 284   for (i = 0; i < N; ++i)
 285     {
 286       gsl_rng_get (r);  /* throw away N iterations */
 287     }
 288
 289   gsl_rng_memcpy (r_save, r);   /* save the intermediate state */
 290
 291   for (i = 0; i < N; ++i)
 292     {
 293       test_a[i] = gsl_rng_get (r);
 294     }
 295
 296   gsl_rng_memcpy (r, r_save);   /* restore the intermediate state */
 297   gsl_rng_free (r_save);
 298
 299   for (i = 0; i < N; ++i)
 300     {
 301       test_b[i] = gsl_rng_get (r);
 302     }
 303
 304   {
 305     int status = 0;
 306     for (i = 0; i < N; ++i)
 307       {
 308         status |= (test_b[i] != test_a[i]);
 309       }
 310     gsl_test (status, "%s, random number state consistency",
 311               gsl_rng_name (r));
 312   }
 313
 314   gsl_rng_free (r);
 315 }
 316
 317
 318 void
 319 rng_parallel_state_test (const gsl_rng_type * T)
 320 {
 321   unsigned long int test_a[N], test_b[N];
 322   unsigned long int test_c[N], test_d[N];
 323   double test_e[N], test_f[N];
 324
 325   int i;
 326
 327   gsl_rng *r1 = gsl_rng_alloc (T);
 328   gsl_rng *r2 = gsl_rng_alloc (T);
 329
 330   for (i = 0; i < N; ++i)
 331     {
 332       gsl_rng_get (r1);         /* throw away N iterations */
 333     }
 334
 335   gsl_rng_memcpy (r2, r1);              /* save the intermediate state */
 336
 337   for (i = 0; i < N; ++i)
 338     {
 339       /* check that there is no hidden state intermixed between r1 and r2 */
 340       test_a[i] = gsl_rng_get (r1);
 341       test_b[i] = gsl_rng_get (r2);
 342       test_c[i] = gsl_rng_uniform_int (r1, 1234);
 343       test_d[i] = gsl_rng_uniform_int (r2, 1234);
 344       test_e[i] = gsl_rng_uniform (r1);
 345       test_f[i] = gsl_rng_uniform (r2);
 346     }
 347
 348   {
 349     int status = 0;
 350     for (i = 0; i < N; ++i)
 351       {
 352         status |= (test_b[i] != test_a[i]);
 353         status |= (test_c[i] != test_d[i]);
 354         status |= (test_e[i] != test_f[i]);
 355       }
 356     gsl_test (status, "%s, parallel random number state consistency",
 357               gsl_rng_name (r1));
 358   }
 359
 360   gsl_rng_free (r1);
 361   gsl_rng_free (r2);
 362
 363 }
 364
 365 void
 366 rng_read_write_test (const gsl_rng_type * T)
 367 {
 368   unsigned long int test_a[N], test_b[N];
 369
 370   int i;
 371
 372   gsl_rng *r = gsl_rng_alloc (T);
 373
 374   for (i = 0; i < N; ++i)
 375     {
 376       gsl_rng_get (r);  /* throw away N iterations */
 377     }
 378
 379   { /* save the state to a binary file */
 380     FILE *f = fopen("test.dat", "wb");
 381     gsl_rng_fwrite(f, r);
 382     fclose(f);
 383   }
 384
 385   for (i = 0; i < N; ++i)
 386     {
 387       test_a[i] = gsl_rng_get (r);
 388     }
 389
 390   { /* read the state from a binary file */
 391     FILE *f = fopen("test.dat", "rb");
 392     gsl_rng_fread(f, r);
 393     fclose(f);
 394   }
 395
 396   for (i = 0; i < N; ++i)
 397     {
 398       test_b[i] = gsl_rng_get (r);
 399     }
 400
 401   {
 402     int status = 0;
 403     for (i = 0; i < N; ++i)
 404       {
 405         status |= (test_b[i] != test_a[i]);
 406       }
 407     gsl_test (status, "%s, random number generator read and write",
 408               gsl_rng_name (r));
 409   }
 410
 411   gsl_rng_free (r);
 412 }
 413
 414 void
 415 generic_rng_test (const gsl_rng_type * T)
 416 {
 417   gsl_rng *r = gsl_rng_alloc (T);
 418   const char *name = gsl_rng_name (r);
 419   unsigned long int kmax = 0, kmin = 1000;
 420   double sigma = 0;
 421   const unsigned long int ran_max = gsl_rng_max (r);
 422   const unsigned long int ran_min = gsl_rng_min (r);
 423
 424   int status = rng_max_test (r, &kmax, ran_max);
 425
 426   gsl_test (status,
 427             "%s, observed vs theoretical maximum (%lu vs %lu)",
 428             name, kmax, ran_max);
 429
 430   status = rng_min_test (r, &kmin, ran_min, ran_max);
 431
 432   gsl_test (status,
 433             "%s, observed vs theoretical minimum (%lu vs %lu)",
 434             name, kmin, ran_min);
 435
 436   status = rng_sum_test (r, &sigma);
 437
 438   gsl_test (status,
 439             "%s, sum test within acceptable sigma (observed %.2g sigma)",
 440             name, sigma);
 441
 442   status = rng_bin_test (r, &sigma);
 443
 444   gsl_test (status,
 445             "%s, bin test within acceptable chisq (observed %.2g sigma)",
 446             name, sigma);
 447
 448   gsl_rng_set (r, 1);   /* set seed to 1 */
 449   status = rng_max_test (r, &kmax, ran_max);
 450
 451   gsl_rng_set (r, 1);   /* set seed to 1 */
 452   status |= rng_min_test (r, &kmin, ran_min, ran_max);
 453
 454   gsl_rng_set (r, 1);   /* set seed to 1 */
 455   status |= rng_sum_test (r, &sigma);
 456
 457   gsl_test (status, "%s, maximum and sum tests for seed=1", name);
 458
 459   gsl_rng_set (r, 12345);       /* set seed to a "typical" value */
 460   status = rng_max_test (r, &kmax, ran_max);
 461
 462   gsl_rng_set (r, 12345);       /* set seed to a "typical" value */
 463   status |= rng_min_test (r, &kmin, ran_min, ran_max);
 464
 465   gsl_rng_set (r, 12345);       /* set seed to a "typical" value */
 466   status |= rng_sum_test (r, &sigma);
 467
 468   gsl_test (status, "%s, maximum and sum tests for non-default seeds", name);
 469
 470   gsl_rng_free (r);
 471 }
 472
 473 int
 474 rng_max_test (gsl_rng * r, unsigned long int *kmax, unsigned long int ran_max)
 475 {
 476   unsigned long int actual_uncovered;
 477   double expect_uncovered;
 478   int status;
 479   unsigned long int max = 0;
 480   int i;
 481
 482   for (i = 0; i < N2; ++i)
 483     {
 484       unsigned long int k = gsl_rng_get (r);
 485       if (k > max)
 486         max = k;
 487     }
 488
 489   *kmax = max;
 490
 491   actual_uncovered = ran_max - max;
 492   expect_uncovered = (double) ran_max / (double) N2;
 493
 494   status = (max > ran_max) || (actual_uncovered > 7 * expect_uncovered) ;
 495
 496   return status;
 497 }
 498
 499 int
 500 rng_min_test (gsl_rng * r, unsigned long int *kmin,
 501               unsigned long int ran_min, unsigned long int ran_max)
 502 {
 503   unsigned long int actual_uncovered;
 504   double expect_uncovered;
 505   int status;
 506   unsigned long int min = 1000000000UL;
 507   int i;
 508
 509   for (i = 0; i < N2; ++i)
 510     {
 511       unsigned long int k = gsl_rng_get (r);
 512       if (k < min)
 513         min = k;
 514     }
 515
 516   *kmin = min;
 517
 518   actual_uncovered = min - ran_min;
 519   expect_uncovered = (double) ran_max / (double) N2;
 520
 521   status = (min < ran_min) || (actual_uncovered > 7 * expect_uncovered);
 522
 523   return status;
 524 }
 525
 526 int
 527 rng_sum_test (gsl_rng * r, double *sigma)
 528 {
 529   double sum = 0;
 530   int i, status;
 531
 532   for (i = 0; i < N2; ++i)
 533     {
 534       double x = gsl_rng_uniform (r) - 0.5;
 535       sum += x;
 536     }
 537
 538   sum /= N2;
 539
 540   /* expect the average to have a variance of 1/(12 n) */
 541
 542   *sigma = sum * sqrt (12.0 * N2);
 543
 544   /* more than 3 sigma is an error (increased to 3.1 to avoid false alarms) */
 545
 546   status = (fabs (*sigma) > 3.1 || fabs(*sigma) < 0.003);
 547
 548   if (status) {
 549       fprintf(stderr,"sum=%g, sigma=%g\n",sum,*sigma);
 550   }
 551
 552   return status;
 553 }
 554
 555 #define BINS 17
 556 #define EXTRA 10
 557 int
 558 rng_bin_test (gsl_rng * r, double *sigma)
 559 {
 560   int count[BINS+EXTRA];
 561   double chisq = 0;
 562   int i, status;
 563
 564   for (i = 0; i < BINS+EXTRA; i++)
 565       count[i] = 0 ;
 566
 567
 568   for (i = 0; i < N2; i++)
 569     {
 570       int j = gsl_rng_uniform_int (r, BINS);
 571       count[j]++ ;
 572     }
 573
 574   chisq = 0 ;
 575   for (i = 0; i < BINS; i++)
 576     {
 577       double x = (double)N2/(double)BINS ;
 578       double d = (count[i] - x) ;
 579       chisq += (d*d) / x;
 580     }
 581
 582   *sigma = sqrt(chisq/BINS) ;
 583
 584   /* more than 3 sigma is an error */
 585
 586   status = (fabs (*sigma) > 3 || fabs(*sigma) < 0.003);
 587
 588   for (i = BINS; i < BINS+EXTRA; i++)
 589     {
 590       if (count[i] != 0)
 591         {
 592           status = 1 ;
 593           gsl_test (status,
 594                     "%s, wrote outside range in bin test "
 595                     "(%d observed vs %d expected)",
 596                     gsl_rng_name(r), i, BINS - 1);
 597         }
 598     }
 599
 600   return status;
 601 }
 602