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

   1 /* rng/random.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 <gsl/gsl_rng.h>
  23
  24 /* This file provides support for random() generators. There are three
  25    versions in widespread use today,
  26
  27    - The original BSD version, e.g. on SunOS 4.1 and FreeBSD.
  28
  29    - The Linux libc5 version, which is differs from the BSD version in
  30    its seeding procedure, possibly due to the introduction of a typo
  31    in the multiplier.
  32
  33    - The GNU glibc2 version, which has a new (and better) seeding
  34    procedure.
  35
  36    They all produce different numbers, due to the different seeding
  37    algorithms, but the algorithm for the generator is the same in each
  38    case.
  39
  40  */
  41
  42 static inline long int random_get (int * i, int * j, int n, long int * x);
  43
  44 static inline unsigned long int random8_get (void *vstate);
  45 static inline unsigned long int random32_get (void *vstate);
  46 static inline unsigned long int random64_get (void *vstate);
  47 static inline unsigned long int random128_get (void *vstate);
  48 static inline unsigned long int random256_get (void *vstate);
  49
  50 static double random8_get_double (void *vstate);
  51 static double random32_get_double (void *vstate);
  52 static double random64_get_double (void *vstate);
  53 static double random128_get_double (void *vstate);
  54 static double random256_get_double (void *vstate);
  55
  56 static void random8_glibc2_set (void *state, unsigned long int s);
  57 static void random32_glibc2_set (void *state, unsigned long int s);
  58 static void random64_glibc2_set (void *state, unsigned long int s);
  59 static void random128_glibc2_set (void *state, unsigned long int s);
  60 static void random256_glibc2_set (void *state, unsigned long int s);
  61
  62 static void random8_libc5_set (void *state, unsigned long int s);
  63 static void random32_libc5_set (void *state, unsigned long int s);
  64 static void random64_libc5_set (void *state, unsigned long int s);
  65 static void random128_libc5_set (void *state, unsigned long int s);
  66 static void random256_libc5_set (void *state, unsigned long int s);
  67
  68 static void random8_bsd_set (void *state, unsigned long int s);
  69 static void random32_bsd_set (void *state, unsigned long int s);
  70 static void random64_bsd_set (void *state, unsigned long int s);
  71 static void random128_bsd_set (void *state, unsigned long int s);
  72 static void random256_bsd_set (void *state, unsigned long int s);
  73
  74 static void bsd_initialize (long int * x, int n, unsigned long int s);
  75 static void libc5_initialize (long int * x, int n, unsigned long int s);
  76 static void glibc2_initialize (long int * x, int n, unsigned long int s);
  77
  78 typedef struct
  79   {
  80     long int x;
  81   }
  82 random8_state_t;
  83
  84 typedef struct
  85   {
  86     int i, j;
  87     long int x[7];
  88   }
  89 random32_state_t;
  90
  91 typedef struct
  92   {
  93     int i, j;
  94     long int x[15];
  95   }
  96 random64_state_t;
  97
  98 typedef struct
  99   {
 100     int i, j;
 101     long int x[31];
 102   }
 103 random128_state_t;
 104
 105 typedef struct
 106   {
 107     int i, j;
 108     long int x[63];
 109   }
 110 random256_state_t;
 111
 112 static inline unsigned long int
 113 random8_get (void *vstate)
 114 {
 115   random8_state_t *state = (random8_state_t *) vstate;
 116
 117   state->x = (1103515245 * state->x + 12345) & 0x7fffffffUL;
 118   return state->x;
 119 }
 120
 121 static inline long int
 122 random_get (int * i, int * j, int n, long int * x)
 123 {
 124   long int k ;
 125
 126   x[*i] += x[*j] ;
 127   k = (x[*i] >> 1) & 0x7FFFFFFF ;
 128
 129   (*i)++ ;
 130   if (*i == n)
 131     *i = 0 ;
 132
 133   (*j)++ ;
 134   if (*j == n)
 135     *j = 0 ;
 136
 137   return k ;
 138 }
 139
 140 static inline unsigned long int
 141 random32_get (void *vstate)
 142 {
 143   random32_state_t *state = (random32_state_t *) vstate;
 144   unsigned long int k = random_get (&state->i, &state->j, 7, state->x) ;
 145   return k ;
 146 }
 147
 148 static inline unsigned long int
 149 random64_get (void *vstate)
 150 {
 151   random64_state_t *state = (random64_state_t *) vstate;
 152   long int k = random_get (&state->i, &state->j, 15, state->x) ;
 153   return k ;
 154 }
 155
 156 static inline unsigned long int
 157 random128_get (void *vstate)
 158 {
 159   random128_state_t *state = (random128_state_t *) vstate;
 160   unsigned long int k = random_get (&state->i, &state->j, 31, state->x) ;
 161   return k ;
 162 }
 163
 164 static inline unsigned long int
 165 random256_get (void *vstate)
 166 {
 167   random256_state_t *state = (random256_state_t *) vstate;
 168   long int k = random_get (&state->i, &state->j, 63, state->x) ;
 169   return k ;
 170 }
 171
 172 static double
 173 random8_get_double (void *vstate)
 174 {
 175   return random8_get (vstate) / 2147483648.0 ;
 176 }
 177
 178 static double
 179 random32_get_double (void *vstate)
 180 {
 181   return random32_get (vstate) / 2147483648.0 ;
 182 }
 183
 184 static double
 185 random64_get_double (void *vstate)
 186 {
 187   return random64_get (vstate) / 2147483648.0 ;
 188 }
 189
 190 static double
 191 random128_get_double (void *vstate)
 192 {
 193   return random128_get (vstate) / 2147483648.0 ;
 194 }
 195
 196 static double
 197 random256_get_double (void *vstate)
 198 {
 199   return random256_get (vstate) / 2147483648.0 ;
 200 }
 201
 202 static void
 203 random8_bsd_set (void *vstate, unsigned long int s)
 204 {
 205   random8_state_t *state = (random8_state_t *) vstate;
 206
 207   if (s == 0)
 208     s = 1;
 209
 210   state->x = s;
 211 }
 212
 213 static void
 214 random32_bsd_set (void *vstate, unsigned long int s)
 215 {
 216   random32_state_t *state = (random32_state_t *) vstate;
 217   int i;
 218
 219   bsd_initialize (state->x, 7, s) ;
 220
 221   state->i = 3;
 222   state->j = 0;
 223
 224   for (i = 0 ; i < 10 * 7 ; i++)
 225     random32_get (state) ;
 226 }
 227
 228 static void
 229 random64_bsd_set (void *vstate, unsigned long int s)
 230 {
 231   random64_state_t *state = (random64_state_t *) vstate;
 232   int i;
 233
 234   bsd_initialize (state->x, 15, s) ;
 235
 236   state->i = 1;
 237   state->j = 0;
 238
 239   for (i = 0 ; i < 10 * 15 ; i++)
 240     random64_get (state) ;
 241 }
 242
 243 static void
 244 random128_bsd_set (void *vstate, unsigned long int s)
 245 {
 246   random128_state_t *state = (random128_state_t *) vstate;
 247   int i;
 248
 249   bsd_initialize (state->x, 31, s) ;
 250
 251   state->i = 3;
 252   state->j = 0;
 253
 254   for (i = 0 ; i < 10 * 31 ; i++)
 255     random128_get (state) ;
 256 }
 257
 258 static void
 259 random256_bsd_set (void *vstate, unsigned long int s)
 260 {
 261   random256_state_t *state = (random256_state_t *) vstate;
 262   int i;
 263
 264   bsd_initialize (state->x, 63, s) ;
 265
 266   state->i = 1;
 267   state->j = 0;
 268
 269   for (i = 0 ; i < 10 * 63 ; i++)
 270     random256_get (state) ;
 271 }
 272
 273 static void
 274 bsd_initialize (long int * x, int n, unsigned long int s)
 275 {
 276   int i;
 277
 278   if (s == 0)
 279     s = 1 ;
 280
 281   x[0] = s;
 282
 283   for (i = 1 ; i < n ; i++)
 284     x[i] = 1103515245 * x[i-1] + 12345 ;
 285 }
 286
 287 static void
 288 libc5_initialize (long int * x, int n, unsigned long int s)
 289 {
 290   int i;
 291
 292   if (s == 0)
 293     s = 1 ;
 294
 295   x[0] = s;
 296
 297   for (i = 1 ; i < n ; i++)
 298     x[i] = 1103515145 * x[i-1] + 12345 ;
 299 }
 300
 301 static void
 302 glibc2_initialize (long int * x, int n, unsigned long int s)
 303 {
 304   int i;
 305
 306   if (s == 0)
 307     s = 1 ;
 308
 309   x[0] = s;
 310
 311   for (i = 1 ; i < n ; i++)
 312     {
 313       const long int h = s / 127773;
 314       const long int t = 16807 * (s - h * 127773) - h * 2836;
 315       if (t < 0)
 316         {
 317           s = t + 2147483647 ;
 318         }
 319       else
 320         {
 321           s = t ;
 322         }
 323
 324     x[i] = s ;
 325     }
 326 }
 327
 328 static void
 329 random8_glibc2_set (void *vstate, unsigned long int s)
 330 {
 331   random8_state_t *state = (random8_state_t *) vstate;
 332
 333   if (s == 0)
 334     s = 1;
 335
 336   state->x = s;
 337 }
 338
 339 static void
 340 random32_glibc2_set (void *vstate, unsigned long int s)
 341 {
 342   random32_state_t *state = (random32_state_t *) vstate;
 343   int i;
 344
 345   glibc2_initialize (state->x, 7, s) ;
 346
 347   state->i = 3;
 348   state->j = 0;
 349
 350   for (i = 0 ; i < 10 * 7 ; i++)
 351     random32_get (state) ;
 352 }
 353
 354 static void
 355 random64_glibc2_set (void *vstate, unsigned long int s)
 356 {
 357   random64_state_t *state = (random64_state_t *) vstate;
 358   int i;
 359
 360   glibc2_initialize (state->x, 15, s) ;
 361
 362   state->i = 1;
 363   state->j = 0;
 364
 365   for (i = 0 ; i < 10 * 15 ; i++)
 366     random64_get (state) ;
 367 }
 368
 369 static void
 370 random128_glibc2_set (void *vstate, unsigned long int s)
 371 {
 372   random128_state_t *state = (random128_state_t *) vstate;
 373   int i;
 374
 375   glibc2_initialize (state->x, 31, s) ;
 376
 377   state->i = 3;
 378   state->j = 0;
 379
 380   for (i = 0 ; i < 10 * 31 ; i++)
 381     random128_get (state) ;
 382 }
 383
 384 static void
 385 random256_glibc2_set (void *vstate, unsigned long int s)
 386 {
 387   random256_state_t *state = (random256_state_t *) vstate;
 388   int i;
 389
 390   glibc2_initialize (state->x, 63, s) ;
 391
 392   state->i = 1;
 393   state->j = 0;
 394
 395   for (i = 0 ; i < 10 * 63 ; i++)
 396     random256_get (state) ;
 397 }
 398
 399
 400 static void
 401 random8_libc5_set (void *vstate, unsigned long int s)
 402 {
 403   random8_state_t *state = (random8_state_t *) vstate;
 404
 405   if (s == 0)
 406     s = 1;
 407
 408   state->x = s;
 409 }
 410
 411 static void
 412 random32_libc5_set (void *vstate, unsigned long int s)
 413 {
 414   random32_state_t *state = (random32_state_t *) vstate;
 415   int i;
 416
 417   libc5_initialize (state->x, 7, s) ;
 418
 419   state->i = 3;
 420   state->j = 0;
 421
 422   for (i = 0 ; i < 10 * 7 ; i++)
 423     random32_get (state) ;
 424 }
 425
 426 static void
 427 random64_libc5_set (void *vstate, unsigned long int s)
 428 {
 429   random64_state_t *state = (random64_state_t *) vstate;
 430   int i;
 431
 432   libc5_initialize (state->x, 15, s) ;
 433
 434   state->i = 1;
 435   state->j = 0;
 436
 437   for (i = 0 ; i < 10 * 15 ; i++)
 438     random64_get (state) ;
 439 }
 440
 441 static void
 442 random128_libc5_set (void *vstate, unsigned long int s)
 443 {
 444   random128_state_t *state = (random128_state_t *) vstate;
 445   int i;
 446
 447   libc5_initialize (state->x, 31, s) ;
 448
 449   state->i = 3;
 450   state->j = 0;
 451
 452   for (i = 0 ; i < 10 * 31 ; i++)
 453     random128_get (state) ;
 454 }
 455
 456 static void
 457 random256_libc5_set (void *vstate, unsigned long int s)
 458 {
 459   random256_state_t *state = (random256_state_t *) vstate;
 460   int i;
 461
 462   libc5_initialize (state->x, 63, s) ;
 463
 464   state->i = 1;
 465   state->j = 0;
 466
 467   for (i = 0 ; i < 10 * 63 ; i++)
 468     random256_get (state) ;
 469 }
 470
 471 static const gsl_rng_type random_glibc2_type =
 472 {"random-glibc2",                       /* name */
 473  0x7fffffffUL,                  /* RAND_MAX */
 474  0,                             /* RAND_MIN */
 475  sizeof (random128_state_t),
 476  &random128_glibc2_set,
 477  &random128_get,
 478  &random128_get_double};
 479
 480 static const gsl_rng_type random8_glibc2_type =
 481 {"random8-glibc2",                      /* name */
 482  0x7fffffffUL,                  /* RAND_MAX */
 483  0,                             /* RAND_MIN */
 484  sizeof (random8_state_t),
 485  &random8_glibc2_set,
 486  &random8_get,
 487  &random8_get_double};
 488
 489 static const gsl_rng_type random32_glibc2_type =
 490 {"random32-glibc2",                     /* name */
 491  0x7fffffffUL,                  /* RAND_MAX */
 492  0,                             /* RAND_MIN */
 493  sizeof (random32_state_t),
 494  &random32_glibc2_set,
 495  &random32_get,
 496  &random32_get_double};
 497
 498 static const gsl_rng_type random64_glibc2_type =
 499 {"random64-glibc2",                     /* name */
 500  0x7fffffffUL,                  /* RAND_MAX */
 501  0,                             /* RAND_MIN */
 502  sizeof (random64_state_t),
 503  &random64_glibc2_set,
 504  &random64_get,
 505  &random64_get_double};
 506
 507 static const gsl_rng_type random128_glibc2_type =
 508 {"random128-glibc2",                    /* name */
 509  0x7fffffffUL,                  /* RAND_MAX */
 510  0,                             /* RAND_MIN */
 511  sizeof (random128_state_t),
 512  &random128_glibc2_set,
 513  &random128_get,
 514  &random128_get_double};
 515
 516 static const gsl_rng_type random256_glibc2_type =
 517 {"random256-glibc2",                    /* name */
 518  0x7fffffffUL,                  /* RAND_MAX */
 519  0,                             /* RAND_MIN */
 520  sizeof (random256_state_t),
 521  &random256_glibc2_set,
 522  &random256_get,
 523  &random256_get_double};
 524
 525 static const gsl_rng_type random_libc5_type =
 526 {"random-libc5",                        /* name */
 527  0x7fffffffUL,                  /* RAND_MAX */
 528  0,                             /* RAND_MIN */
 529  sizeof (random128_state_t),
 530  &random128_libc5_set,
 531  &random128_get,
 532  &random128_get_double};
 533
 534 static const gsl_rng_type random8_libc5_type =
 535 {"random8-libc5",                       /* name */
 536  0x7fffffffUL,                  /* RAND_MAX */
 537  0,                             /* RAND_MIN */
 538  sizeof (random8_state_t),
 539  &random8_libc5_set,
 540  &random8_get,
 541  &random8_get_double};
 542
 543 static const gsl_rng_type random32_libc5_type =
 544 {"random32-libc5",                      /* name */
 545  0x7fffffffUL,                  /* RAND_MAX */
 546  0,                             /* RAND_MIN */
 547  sizeof (random32_state_t),
 548  &random32_libc5_set,
 549  &random32_get,
 550  &random32_get_double};
 551
 552 static const gsl_rng_type random64_libc5_type =
 553 {"random64-libc5",                      /* name */
 554  0x7fffffffUL,                  /* RAND_MAX */
 555  0,                             /* RAND_MIN */
 556  sizeof (random64_state_t),
 557  &random64_libc5_set,
 558  &random64_get,
 559  &random64_get_double};
 560
 561 static const gsl_rng_type random128_libc5_type =
 562 {"random128-libc5",                     /* name */
 563  0x7fffffffUL,                  /* RAND_MAX */
 564  0,                             /* RAND_MIN */
 565  sizeof (random128_state_t),
 566  &random128_libc5_set,
 567  &random128_get,
 568  &random128_get_double};
 569
 570 static const gsl_rng_type random256_libc5_type =
 571 {"random256-libc5",                     /* name */
 572  0x7fffffffUL,                  /* RAND_MAX */
 573  0,                             /* RAND_MIN */
 574  sizeof (random256_state_t),
 575  &random256_libc5_set,
 576  &random256_get,
 577  &random256_get_double};
 578
 579 static const gsl_rng_type random_bsd_type =
 580 {"random-bsd",                  /* name */
 581  0x7fffffffUL,                  /* RAND_MAX */
 582  0,                             /* RAND_MIN */
 583  sizeof (random128_state_t),
 584  &random128_bsd_set,
 585  &random128_get,
 586  &random128_get_double};
 587
 588 static const gsl_rng_type random8_bsd_type =
 589 {"random8-bsd",                 /* name */
 590  0x7fffffffUL,                  /* RAND_MAX */
 591  0,                             /* RAND_MIN */
 592  sizeof (random8_state_t),
 593  &random8_bsd_set,
 594  &random8_get,
 595  &random8_get_double};
 596
 597 static const gsl_rng_type random32_bsd_type =
 598 {"random32-bsd",                        /* name */
 599  0x7fffffffUL,                  /* RAND_MAX */
 600  0,                             /* RAND_MIN */
 601  sizeof (random32_state_t),
 602  &random32_bsd_set,
 603  &random32_get,
 604  &random32_get_double};
 605
 606 static const gsl_rng_type random64_bsd_type =
 607 {"random64-bsd",                        /* name */
 608  0x7fffffffUL,                  /* RAND_MAX */
 609  0,                             /* RAND_MIN */
 610  sizeof (random64_state_t),
 611  &random64_bsd_set,
 612  &random64_get,
 613  &random64_get_double};
 614
 615 static const gsl_rng_type random128_bsd_type =
 616 {"random128-bsd",               /* name */
 617  0x7fffffffUL,                  /* RAND_MAX */
 618  0,                             /* RAND_MIN */
 619  sizeof (random128_state_t),
 620  &random128_bsd_set,
 621  &random128_get,
 622  &random128_get_double};
 623
 624 static const gsl_rng_type random256_bsd_type =
 625 {"random256-bsd",               /* name */
 626  0x7fffffffUL,                  /* RAND_MAX */
 627  0,                             /* RAND_MIN */
 628  sizeof (random256_state_t),
 629  &random256_bsd_set,
 630  &random256_get,
 631  &random256_get_double};
 632
 633 const gsl_rng_type *gsl_rng_random_libc5    = &random_libc5_type;
 634 const gsl_rng_type *gsl_rng_random8_libc5   = &random8_libc5_type;
 635 const gsl_rng_type *gsl_rng_random32_libc5  = &random32_libc5_type;
 636 const gsl_rng_type *gsl_rng_random64_libc5  = &random64_libc5_type;
 637 const gsl_rng_type *gsl_rng_random128_libc5 = &random128_libc5_type;
 638 const gsl_rng_type *gsl_rng_random256_libc5 = &random256_libc5_type;
 639
 640 const gsl_rng_type *gsl_rng_random_glibc2    = &random_glibc2_type;
 641 const gsl_rng_type *gsl_rng_random8_glibc2   = &random8_glibc2_type;
 642 const gsl_rng_type *gsl_rng_random32_glibc2  = &random32_glibc2_type;
 643 const gsl_rng_type *gsl_rng_random64_glibc2  = &random64_glibc2_type;
 644 const gsl_rng_type *gsl_rng_random128_glibc2 = &random128_glibc2_type;
 645 const gsl_rng_type *gsl_rng_random256_glibc2 = &random256_glibc2_type;
 646
 647 const gsl_rng_type *gsl_rng_random_bsd    = &random_bsd_type;
 648 const gsl_rng_type *gsl_rng_random8_bsd   = &random8_bsd_type;
 649 const gsl_rng_type *gsl_rng_random32_bsd  = &random32_bsd_type;
 650 const gsl_rng_type *gsl_rng_random64_bsd  = &random64_bsd_type;
 651 const gsl_rng_type *gsl_rng_random128_bsd = &random128_bsd_type;
 652 const gsl_rng_type *gsl_rng_random256_bsd = &random256_bsd_type;
 653
 654
 655
 656