git.cworth.org Git - vogl/blob - src/voglcore/vogl_rand.cpp

   1 /**************************************************************************
   2  *
   3  * Copyright 2013-2014 RAD Game Tools and Valve Software
   4  * Copyright 2010-2014 Rich Geldreich and Tenacious Software LLC
   5  * All Rights Reserved.
   6  *
   7  * Permission is hereby granted, free of charge, to any person obtaining a copy
   8  * of this software and associated documentation files (the "Software"), to deal
   9  * in the Software without restriction, including without limitation the rights
  10  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  11  * copies of the Software, and to permit persons to whom the Software is
  12  * furnished to do so, subject to the following conditions:
  13  *
  14  * The above copyright notice and this permission notice shall be included in
  15  * all copies or substantial portions of the Software.
  16  *
  17  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  18  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  19  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  20  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  21  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  22  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  23  * THE SOFTWARE.
  24  *
  25  **************************************************************************/
  26
  27 // File: vogl_rand.cpp
  28 // See:
  29 // http://www.ciphersbyritter.com/NEWS4/RANDC.HTM
  30 // http://burtleburtle.net/bob/rand/smallprng.html
  31 // http://www.cs.ucl.ac.uk/staff/d.jones/GoodPracticeRNG.pdf
  32 // See GPG7, page 120, or http://www.lomont.org/Math/Papers/2008/Lomont_PRNG_2008.pdf
  33 // I've tested random::urand32() using the CISC test suite here: http://www.csis.hku.hk/cisc/download/idetect/
  34 #include "vogl_core.h"
  35 #include "vogl_rand.h"
  36 #include "vogl_hash.h"
  37 #include "vogl_md5.h"
  38 #include "vogl_bigint128.h"
  39
  40 #include "vogl_image_utils.h"
  41
  42 // TODO: Linux specific, needed by the seed_from_urand() method
  43 #include <unistd.h>
  44 #include <sys/syscall.h>
  45
  46 #define znew (z = 36969 * (z & 65535) + (z >> 16))
  47 #define wnew (w = 18000 * (w & 65535) + (w >> 16))
  48 #define MWC ((znew << 16) + wnew)
  49 #define FIB ((b = a + b), (a = b - a))
  50 #define KISS ((MWC ^ VOGL_CONG) + VOGL_SHR3)
  51 #define LFIB4 (c++, t[c] = t[c] + t[UC(c + 58)] + t[UC(c + 119)] + t[UC(c + 178)])
  52 #define SWB (c++, bro = (x < y), t[c] = (x = t[UC(c + 34)]) - (y = t[UC(c + 19)] + bro))
  53 #define UNI (KISS * 2.328306e-10)
  54 #define VNI ((long)KISS) * 4.656613e-10
  55 #define UC (unsigned char) /*a cast operation*/
  56
  57 //#define ROTL(x,k) (((x)<<(k))|((x)>>(32-(k))))
  58 #define ROTL(x, k) VOGL_ROTATE_LEFT(x, k)
  59
  60 namespace vogl
  61 {
  62     random g_random;
  63     fast_random g_fast_random;
  64     thread_safe_random g_thread_safe_random;
  65
  66     static const double cNorm32 = 1.0 / (double)0x100000000ULL;
  67
  68     kiss99::kiss99()
  69     {
  70         x = 123456789;
  71         y = 362436000;
  72         z = 521288629;
  73         c = 7654321;
  74     }
  75
  76     void kiss99::seed(uint32 i, uint32 j, uint32 k)
  77     {
  78         x = i;
  79         y = j;
  80         z = k;
  81         c = 7654321;
  82     }
  83
  84     inline uint32 kiss99::next()
  85     {
  86         x = 69069 * x + 12345;
  87
  88         y ^= (y << 13);
  89         y ^= (y >> 17);
  90         y ^= (y << 5);
  91
  92         uint64_t t = c;
  93         t += (698769069ULL * z);
  94         c = static_cast<uint32>(t >> 32);
  95         z = static_cast<uint32>(t);
  96
  97         return (x + y + z);
  98     }
  99
 100     inline uint32 ranctx::next()
 101     {
 102         uint32 e = a - ROTL(b, 27);
 103         a = b ^ ROTL(c, 17);
 104         b = c + d;
 105         c = d + e;
 106         d = e + a;
 107         return d;
 108     }
 109
 110     void ranctx::seed(uint32 seed)
 111     {
 112         a = 0xf1ea5eed, b = c = d = seed;
 113         for (uint32 i = 0; i < 20; ++i)
 114             next();
 115     }
 116
 117     well512::well512()
 118     {
 119         seed(0xDEADBE3F);
 120     }
 121
 122     void well512::seed(uint32 seed[well512::cStateSize])
 123     {
 124         memcpy(m_state, seed, sizeof(m_state));
 125         m_index = 0;
 126     }
 127
 128     void well512::seed(uint32 seed)
 129     {
 130         uint32 jsr = utils::swap32(seed) ^ 0xAAC29377;
 131
 132         for (uint i = 0; i < cStateSize; i++)
 133         {
 134             VOGL_SHR3;
 135             seed = bitmix32c(seed);
 136
 137             m_state[i] = seed ^ jsr;
 138         }
 139         m_index = 0;
 140     }
 141
 142     void well512::seed(uint32 seed1, uint32 seed2, uint32 seed3)
 143     {
 144         uint32 jsr = seed2;
 145         uint32 jcong = seed3;
 146
 147         for (uint i = 0; i < cStateSize; i++)
 148         {
 149             VOGL_SHR3;
 150             seed1 = bitmix32c(seed1);
 151             VOGL_CONG;
 152
 153             m_state[i] = seed1 ^ jsr ^ jcong;
 154         }
 155         m_index = 0;
 156     }
 157
 158     void well512::seed(uint32 seed1, uint32 seed2, uint32 seed3, uint seed4)
 159     {
 160         uint32 jsr = seed2;
 161         uint32 jcong = seed3;
 162
 163         m_state[0] = seed1;
 164         m_state[1] = seed2;
 165         m_state[2] = seed3;
 166         m_state[3] = seed4;
 167
 168         for (uint i = 4; i < cStateSize; i++)
 169         {
 170             VOGL_SHR3;
 171             VOGL_CONG;
 172
 173             m_state[i] = jsr ^ jcong ^ m_state[(i - 4) & (cStateSize - 1)];
 174         }
 175
 176         m_index = 0;
 177     }
 178
 179     inline uint32 well512::next()
 180     {
 181         uint32 a, b, c, d;
 182         a = m_state[m_index];
 183         c = m_state[(m_index + 13) & 15];
 184         b = a ^ c ^ (a << 16) ^ (c << 15);
 185         c = m_state[(m_index + 9) & 15];
 186         c ^= (c >> 11);
 187         a = m_state[m_index] = b ^ c;
 188         d = a ^ ((a << 5) & 0xDA442D20UL);
 189         m_index = (m_index + 15) & 15;
 190         a = m_state[m_index];
 191         m_state[m_index] = a ^ b ^ d ^ (a << 2) ^ (b << 18) ^ (c << 28);
 192         return m_state[m_index];
 193     }
 194
 195     random::random()
 196     {
 197         seed(12345, 65435, 34221);
 198     }
 199
 200     random::random(uint32 i)
 201     {
 202         seed(i);
 203     }
 204
 205     void random::seed()
 206     {
 207         seed(12345, 65435, 34221);
 208     }
 209
 210     void random::seed(uint32 i1, uint32 i2, uint32 i3)
 211     {
 212         m_ranctx.seed(i1 ^ i2 ^ i3);
 213
 214         m_kiss99.seed(i1, i2, i3);
 215
 216         m_well512.seed(i1, i2, i3);
 217
 218         m_pdes.seed((static_cast<uint64_t>(i1) | (static_cast<uint64_t>(i2) << 32U)) ^ (static_cast<uint64_t>(i3) << 10));
 219
 220         const uint warming_iters = 16;
 221         for (uint i = 0; i < warming_iters; i++)
 222             urand32();
 223     }
 224
 225     void random::seed(uint32 i1, uint32 i2, uint32 i3, uint i4)
 226     {
 227         m_ranctx.seed(i1 ^ i2 ^ i3);
 228
 229         m_kiss99.seed(i1, i2, i4);
 230
 231         m_well512.seed(i1, i2, i3, i4);
 232
 233         m_pdes.seed((static_cast<uint64_t>(i1) | (static_cast<uint64_t>(i4) << 32U)) ^ (static_cast<uint64_t>(i3) << 10));
 234
 235         const uint warming_iters = 16;
 236         for (uint i = 0; i < warming_iters; i++)
 237             urand32();
 238     }
 239
 240     void random::seed(const md5_hash &h)
 241     {
 242         seed(h[0], h[1], h[2], h[3]);
 243     }
 244
 245     void random::seed(uint32 i)
 246     {
 247         md5_hash h(md5_hash_gen(&i, sizeof(i)).finalize());
 248         seed(h[0], h[1], h[2], h[3]);
 249     }
 250
 251     void random::seed_from_urandom()
 252     {
 253         const uint N = 4;
 254         uint32 buf[N] = { 0xABCDEF, 0x12345678, 0xFFFECABC, 0xABCDDEF0 };
 255         FILE *fp = vogl_fopen("/dev/urandom", "rb");
 256         if (fp)
 257         {
 258             size_t n = fread(buf, 1, sizeof(buf), fp);
 259             VOGL_NOTE_UNUSED(n);
 260             vogl_fclose(fp);
 261         }
 262
 263         // Hash thread, process ID's, etc. for good measure, or in case the call failed, or to make this easier to port. (we can just leave out the /dev/urandom read)
 264         pid_t tid = (pid_t)syscall(SYS_gettid);
 265         pid_t pid = getpid();
 266         pid_t ppid = getppid();
 267
 268         buf[0] ^= static_cast<uint32>(time(NULL));
 269         buf[1] ^= static_cast<uint32>(utils::RDTSC()) ^ static_cast<uint32>(ppid);
 270         buf[2] ^= static_cast<uint32>(tid);
 271         buf[3] ^= static_cast<uint32>(pid);
 272
 273         seed(buf[0], buf[1], buf[2], buf[3]);
 274     }
 275
 276     uint32 random::urand32()
 277     {
 278         return m_kiss99.next() ^ m_ranctx.next() ^ m_well512.next() ^ m_pdes.next32();
 279     }
 280
 281     uint64_t random::urand64()
 282     {
 283         uint64_t result = urand32();
 284         result <<= 32ULL;
 285         result |= urand32();
 286         return result;
 287     }
 288
 289     uint32 random::fast_urand32()
 290     {
 291         return m_well512.next();
 292     }
 293
 294     uint32 random::get_bit()
 295     {
 296         uint32 k = urand32();
 297         return (k ^ (k >> 6) ^ (k >> 10) ^ (k >> 30)) & 1;
 298     }
 299
 300     double random::drand(double l, double h)
 301     {
 302         VOGL_ASSERT(l <= h);
 303         if (l >= h)
 304             return l;
 305
 306         double fract = (urand32() * (cNorm32 * cNorm32)) + (urand32() * cNorm32);
 307
 308         // clamp here shouldn't be necessary, but it ensures this function provably cannot return out of range results
 309         return math::clamp(l + (h - l) * fract, l, h);
 310     }
 311
 312     float random::frand(float l, float h)
 313     {
 314         VOGL_ASSERT(l <= h);
 315         if (l >= h)
 316             return l;
 317
 318         float r = static_cast<float>(l + (h - l) * (urand32() * cNorm32));
 319
 320         // clamp here shouldn't be necessary, but it ensures this function provably cannot return out of range results
 321         return math::clamp<float>(r, l, h);
 322     }
 323
 324     static inline uint umul32_return_high(uint32 range, uint32 rnd)
 325     {
 326 #if defined(_M_IX86) && defined(_MSC_VER)
 327         //uint32 rnd_range = static_cast<uint32>(__emulu(range, rnd) >> 32U);
 328         uint32 x[2];
 329         *reinterpret_cast<uint64_t *>(x) = __emulu(range, rnd);
 330         uint32 rnd_range = x[1];
 331 #else
 332         uint32 rnd_range = static_cast<uint32>((range * static_cast<uint64_t>(rnd)) >> 32U);
 333 #endif
 334
 335         return rnd_range;
 336     }
 337
 338     int random::irand(int l, int h)
 339     {
 340         VOGL_ASSERT(l < h);
 341         if (l >= h)
 342             return l;
 343
 344         uint32 range = static_cast<uint32>(h - l);
 345
 346         uint32 rnd = urand32();
 347
 348         uint32 rnd_range = umul32_return_high(range, rnd);
 349
 350         int result = l + rnd_range;
 351         VOGL_ASSERT((result >= l) && (result < h));
 352         return result;
 353     }
 354
 355     int random::irand_inclusive(int l, int h)
 356     {
 357         VOGL_ASSERT(l <= h);
 358         if (l >= h)
 359             return l;
 360
 361         uint32 range = static_cast<uint32>(h - l);
 362
 363         uint32 rnd = urand32();
 364
 365         int result = static_cast<int>(rnd);
 366         if (range != cUINT32_MAX)
 367         {
 368             ++range;
 369             uint32 rnd_range = umul32_return_high(range, rnd);
 370             result = l + rnd_range;
 371         }
 372
 373         VOGL_ASSERT((result >= l) && (result <= h));
 374         return result;
 375     }
 376
 377     uint random::urand(uint l, uint h)
 378     {
 379         VOGL_ASSERT(l < h);
 380         if (l >= h)
 381             return l;
 382
 383         uint32 range = h - l;
 384
 385         uint32 rnd = urand32();
 386
 387         uint32 rnd_range = umul32_return_high(range, rnd);
 388
 389         uint result = l + rnd_range;
 390         VOGL_ASSERT((result >= l) && (result < h));
 391         return result;
 392     }
 393
 394     uint random::urand_inclusive(uint l, uint h)
 395     {
 396         VOGL_ASSERT(l <= h);
 397         if (l >= h)
 398             return l;
 399
 400         uint32 range = h - l;
 401
 402         uint32 rnd = urand32();
 403
 404         uint32 result = rnd;
 405         if (range != cUINT32_MAX)
 406         {
 407             ++range;
 408             uint32 rnd_range = umul32_return_high(range, rnd);
 409             result = l + rnd_range;
 410         }
 411
 412         VOGL_ASSERT((result >= l) && (result <= h));
 413         return result;
 414     }
 415
 416     uint64_t random::urand64(uint64_t l, uint64_t h)
 417     {
 418         VOGL_ASSERT(l < h);
 419         if (l >= h)
 420             return l;
 421
 422         uint64_t range = h - l;
 423
 424         uint64_t rnd = urand64();
 425
 426         bigint128 mul_result;
 427         bigint128::unsigned_multiply(range, rnd, mul_result);
 428
 429         uint64_t rnd_range = mul_result.get_qword(1);
 430
 431         uint64_t result = l + rnd_range;
 432         VOGL_ASSERT((result >= l) && (result < h));
 433         return result;
 434     }
 435
 436     uint64_t random::urand64_inclusive(uint64_t l, uint64_t h)
 437     {
 438         VOGL_ASSERT(l <= h);
 439         if (l >= h)
 440             return l;
 441
 442         uint64_t range = h - l;
 443
 444         uint64_t rnd = urand64();
 445
 446         uint64_t result = rnd;
 447         if (range != cUINT64_MAX)
 448         {
 449             bigint128 mul_result;
 450             bigint128::unsigned_multiply(range + 1, rnd, mul_result);
 451
 452             uint64_t rnd_range = mul_result.get_qword(1);
 453
 454             result = l + rnd_range;
 455         }
 456
 457         VOGL_ASSERT((result >= l) && (result <= h));
 458         return result;
 459     }
 460
 461     /*
 462    ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
 463    THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
 464    VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.
 465    The function returns a normally distributed pseudo-random number
 466    with a given mean and standard devaiation.  Calls are made to a
 467    function subprogram which must return independent random
 468    numbers uniform in the interval (0,1).
 469    The algorithm uses the ratio of uniforms method of A.J. Kinderman
 470    and J.F. Monahan augmented with quadratic bounding curves.
 471    */
 472     double random::gaussian(double mean, double stddev)
 473     {
 474         double q, u, v, x, y;
 475
 476         /*
 477         Generate P = (u,v) uniform in rect. enclosing acceptance region
 478         Make sure that any random numbers <= 0 are rejected, since
 479         gaussian() requires uniforms > 0, but RandomUniform() delivers >= 0.
 480         */
 481         do
 482         {
 483             u = drand(0, 1);
 484             v = drand(0, 1);
 485             if (u <= 0.0 || v <= 0.0)
 486             {
 487                 u = 1.0;
 488                 v = 1.0;
 489             }
 490             v = 1.7156 * (v - 0.5);
 491
 492             /*  Evaluate the quadratic form */
 493             x = u - 0.449871;
 494             y = fabs(v) + 0.386595;
 495             q = x * x + y * (0.19600 * y - 0.25472 * x);
 496
 497             /* Accept P if inside inner ellipse */
 498             if (q < 0.27597)
 499                 break;
 500
 501             /*  Reject P if outside outer ellipse, or outside acceptance region */
 502         } while ((q > 0.27846) || (v * v > -4.0 * log(u) * u * u));
 503
 504         /*  Return ratio of P's coordinates as the normal deviate */
 505         return (mean + stddev * v / u);
 506     }
 507
 508     thread_safe_random::thread_safe_random()
 509     {
 510         scoped_spinlock scoped_lock(m_spinlock);
 511         m_rm.seed();
 512     }
 513
 514     thread_safe_random::thread_safe_random(uint32 i)
 515     {
 516         scoped_spinlock scoped_lock(m_spinlock);
 517         m_rm.seed(i);
 518     }
 519
 520     void thread_safe_random::seed()
 521     {
 522         scoped_spinlock scoped_lock(m_spinlock);
 523         m_rm.seed();
 524     }
 525
 526     void thread_safe_random::seed(uint32 i)
 527     {
 528         scoped_spinlock scoped_lock(m_spinlock);
 529         m_rm.seed(i);
 530     }
 531
 532     void thread_safe_random::seed(uint32 i1, uint32 i2, uint32 i3)
 533     {
 534         scoped_spinlock scoped_lock(m_spinlock);
 535         m_rm.seed(i1, i2, i3);
 536     }
 537
 538     void thread_safe_random::seed(uint32 i1, uint32 i2, uint32 i3, uint32 i4)
 539     {
 540         scoped_spinlock scoped_lock(m_spinlock);
 541         m_rm.seed(i1, i2, i3, i4);
 542     }
 543
 544     void thread_safe_random::seed(const md5_hash &h)
 545     {
 546         scoped_spinlock scoped_lock(m_spinlock);
 547         m_rm.seed(h);
 548     }
 549
 550     void thread_safe_random::seed_from_urandom()
 551     {
 552         scoped_spinlock scoped_lock(m_spinlock);
 553         m_rm.seed_from_urandom();
 554     }
 555
 556     uint8 thread_safe_random::urand8()
 557     {
 558         scoped_spinlock scoped_lock(m_spinlock);
 559         return m_rm.urand8();
 560     }
 561
 562     uint16 thread_safe_random::urand16()
 563     {
 564         scoped_spinlock scoped_lock(m_spinlock);
 565         return m_rm.urand16();
 566     }
 567
 568     uint32 thread_safe_random::urand32()
 569     {
 570         scoped_spinlock scoped_lock(m_spinlock);
 571         return m_rm.urand32();
 572     }
 573
 574     uint64_t thread_safe_random::urand64()
 575     {
 576         scoped_spinlock scoped_lock(m_spinlock);
 577         return m_rm.urand64();
 578     }
 579
 580     uint32 thread_safe_random::fast_urand32()
 581     {
 582         scoped_spinlock scoped_lock(m_spinlock);
 583         return m_rm.fast_urand32();
 584     }
 585
 586     uint32 thread_safe_random::get_bit()
 587     {
 588         scoped_spinlock scoped_lock(m_spinlock);
 589         return m_rm.get_bit();
 590     }
 591
 592     double thread_safe_random::drand(double l, double h)
 593     {
 594         scoped_spinlock scoped_lock(m_spinlock);
 595         return m_rm.drand(l, h);
 596     }
 597
 598     float thread_safe_random::frand(float l, float h)
 599     {
 600         scoped_spinlock scoped_lock(m_spinlock);
 601         return m_rm.frand(l, h);
 602     }
 603
 604     int thread_safe_random::irand(int l, int h)
 605     {
 606         scoped_spinlock scoped_lock(m_spinlock);
 607         return m_rm.irand(l, h);
 608     }
 609
 610     int thread_safe_random::irand_inclusive(int l, int h)
 611     {
 612         scoped_spinlock scoped_lock(m_spinlock);
 613         return m_rm.irand_inclusive(l, h);
 614     }
 615
 616     uint thread_safe_random::urand(uint l, uint h)
 617     {
 618         scoped_spinlock scoped_lock(m_spinlock);
 619         return m_rm.urand(l, h);
 620     }
 621
 622     uint thread_safe_random::urand_inclusive(uint l, uint h)
 623     {
 624         scoped_spinlock scoped_lock(m_spinlock);
 625         return m_rm.urand_inclusive(l, h);
 626     }
 627
 628     // Returns random 64-bit unsigned int between [l, h)
 629     uint64_t thread_safe_random::urand64(uint64_t l, uint64_t h)
 630     {
 631         scoped_spinlock scoped_lock(m_spinlock);
 632         return m_rm.urand64(l, h);
 633     }
 634
 635     // Returns random 64-bit unsigned int between [l, h]
 636     uint64_t thread_safe_random::urand64_inclusive(uint64_t l, uint64_t h)
 637     {
 638         scoped_spinlock scoped_lock(m_spinlock);
 639         return m_rm.urand64_inclusive(l, h);
 640     }
 641
 642     double thread_safe_random::gaussian(double mean, double stddev)
 643     {
 644         scoped_spinlock scoped_lock(m_spinlock);
 645         return m_rm.gaussian(mean, stddev);
 646     }
 647
 648     int fast_random::irand(int l, int h)
 649     {
 650         VOGL_ASSERT(l < h);
 651         if (l >= h)
 652             return l;
 653
 654         uint32 range = static_cast<uint32>(h - l);
 655
 656         uint32 rnd = urand32();
 657
 658         uint32 rnd_range = umul32_return_high(range, rnd);
 659
 660         int result = l + rnd_range;
 661         VOGL_ASSERT((result >= l) && (result < h));
 662         return result;
 663     }
 664
 665     uint fast_random::urand(uint l, uint h)
 666     {
 667         VOGL_ASSERT(l < h);
 668         if (l >= h)
 669             return l;
 670
 671         uint32 range = h - l;
 672
 673         uint32 rnd = urand32();
 674
 675         uint32 rnd_range = umul32_return_high(range, rnd);
 676
 677         uint result = l + rnd_range;
 678         VOGL_ASSERT((result >= l) && (result < h));
 679         return result;
 680     }
 681
 682     double fast_random::drand(double l, double h)
 683     {
 684         VOGL_ASSERT(l <= h);
 685         if (l >= h)
 686             return l;
 687
 688         return math::clamp(l + (h - l) * (urand32() * cNorm32), l, h);
 689     }
 690
 691     float fast_random::frand(float l, float h)
 692     {
 693         VOGL_ASSERT(l <= h);
 694         if (l >= h)
 695             return l;
 696
 697         float r = static_cast<float>(l + (h - l) * (urand32() * cNorm32));
 698
 699         return math::clamp<float>(r, l, h);
 700     }
 701
 702     // Got this from http://www.physics.buffalo.edu/phy410-505-2009/topic2/lec-2-1.pdf
 703     template <typename T>
 704     double chi_sqr_test(T gen)
 705     {
 706         int tries = 10000000;
 707         const int BINS = 10000;
 708         printf("Binning %i tries in %i bins\n", tries, BINS);
 709         int bin[BINS];
 710         for (int b = 0; b < BINS; b++)
 711             bin[b] = 0;
 712
 713         for (int i = 0; i < tries; i++)
 714         {
 715             int b = gen(BINS);
 716             VOGL_ASSERT(b < BINS);
 717             if (b == BINS)
 718                 --b;
 719             ++bin[b];
 720         }
 721         double chiSqr = 0;
 722         double expect = tries / double(BINS);
 723         for (int b = 0; b < BINS; b++)
 724         {
 725             double diff = bin[b] - expect;
 726             chiSqr += diff * diff;
 727         }
 728         chiSqr /= expect;
 729         printf("chi-square: %f, chi-square/d.o.f. = %f\n", chiSqr, chiSqr / (BINS - 1));
 730         return chiSqr / (BINS - 1);
 731     }
 732
 733     struct fast_rand_test_obj
 734     {
 735         fast_rand_test_obj(fast_random &r)
 736             : m_r(r)
 737         {
 738             printf("%s:\n", VOGL_METHOD_NAME);
 739         }
 740         int operator()(int limit)
 741         {
 742             return m_r.irand(0, limit);
 743         }
 744         fast_random &m_r;
 745     };
 746
 747     struct rand_test_obj
 748     {
 749         rand_test_obj(random &r)
 750             : m_r(r)
 751         {
 752             printf("%s:\n", VOGL_METHOD_NAME);
 753         }
 754         int operator()(int limit)
 755         {
 756             return m_r.irand(0, limit);
 757         }
 758         random &m_r;
 759     };
 760
 761     struct rand_test_obj2
 762     {
 763         rand_test_obj2(random &r)
 764             : m_r(r)
 765         {
 766             printf("%s:\n", VOGL_METHOD_NAME);
 767         }
 768         int operator()(int limit)
 769         {
 770             uint b = 0;
 771             for (uint i = 0; i < 32; i++)
 772                 b = (b << 1) | m_r.get_bit();
 773             return (static_cast<double>(b) * limit) / (static_cast<double>(cUINT32_MAX) + 1.0f);
 774         }
 775         random &m_r;
 776     };
 777
 778     struct dbl_rand_test_obj
 779     {
 780         dbl_rand_test_obj(random &r)
 781             : m_r(r)
 782         {
 783             printf("%s:\n", VOGL_METHOD_NAME);
 784         }
 785         int operator()(int limit)
 786         {
 787             return static_cast<int>(m_r.drand(0, 1.0f) * limit);
 788         }
 789         random &m_r;
 790     };
 791
 792     template <typename T>
 793     double pi_estimate(const unsigned long points, T &func)
 794     {
 795         unsigned long hits = 0;
 796         for (unsigned long i = 0; i < points; i++)
 797         {
 798             double x = func.drand(0.0f, 1.0f);
 799             double y = func.drand(0.0f, 1.0f);
 800             if ((x - 0.5) * (x - 0.5) + (y - 0.5) * (y - 0.5) < 0.25)
 801                 ++hits;
 802         }
 803         return 4.0f * double(hits) / double(points);
 804     }
 805
 806     template <typename T>
 807     void measure_pi(
 808         const unsigned long trials,
 809         const unsigned long points,
 810         double &average,
 811         double &std_dev, T &func)
 812     {
 813         double sum = 0;
 814         double squared_sum = 0;
 815         for (unsigned long t = 0; t < trials; t++)
 816         {
 817             double pi = pi_estimate(points, func);
 818             sum += pi;
 819             squared_sum += pi * pi;
 820         }
 821         average = sum / trials;
 822         std_dev = math::maximum<double>(0.0f, squared_sum / trials - average * average);
 823         std_dev = sqrt(std_dev / (trials - 1));
 824     }
 825
 826     // Very basic tests, but better than nothing.
 827     bool rand_test()
 828     {
 829         fast_random frm;
 830         random rm;
 831
 832         time_t seconds = time(NULL);
 833         uint32 seed = static_cast<uint32>(seconds);
 834         frm.seed(seed);
 835         rm.seed(seed);
 836
 837         image_u8 eye_test_img(1024, 1024);
 838         for (uint t = 0; t < 500; t++)
 839         {
 840             for (uint i = 0; i < 40000; i++)
 841             {
 842                 uint x, y;
 843                 if (i > 10000)
 844                 {
 845                     x = static_cast<int>(rm.frand(0, 1.0f) * eye_test_img.get_width());
 846                     y = static_cast<int>(rm.frand(0, 1.0f) * eye_test_img.get_height());
 847                 }
 848                 else
 849                 {
 850                     x = rm.irand(0, eye_test_img.get_width());
 851                     y = rm.irand(0, eye_test_img.get_height());
 852                 }
 853                 uint l = rm.irand(1, 16);
 854                 color_quad_u8 c(eye_test_img(x, y));
 855                 int k = math::minimum<int>(c[0] + l, 255);
 856                 c.set(k);
 857                 eye_test_img(x, y) = c;
 858             }
 859
 860             dynamic_string filename(cVarArg, "rand_eye_test_%02i.png", t);
 861             image_utils::write_to_file(filename.get_ptr(), eye_test_img, image_utils::cWriteFlagIgnoreAlpha);
 862             printf("Wrote %s\n", filename.get_ptr());
 863         }
 864
 865         chi_sqr_test(fast_rand_test_obj(frm));
 866         chi_sqr_test(fast_rand_test_obj(frm));
 867
 868         chi_sqr_test(rand_test_obj(rm));
 869         chi_sqr_test(rand_test_obj(rm));
 870         chi_sqr_test(rand_test_obj2(rm));
 871
 872         chi_sqr_test(dbl_rand_test_obj(rm));
 873         chi_sqr_test(dbl_rand_test_obj(rm));
 874
 875         for (uint i = 0; i < 5; i++)
 876         {
 877             double avg, std_dev;
 878             measure_pi(10000, 10000, avg, std_dev, rm);
 879             printf("pi avg: %1.15f, std dev: %1.15f\n", avg, std_dev);
 880
 881             // check for obviously bogus results
 882             if ((fabs(avg) - 3.14f) > .01f)
 883                 return false;
 884
 885             if (fabs(std_dev) < .0001f)
 886                 return false;
 887         }
 888
 889         random cur_rm(rm);
 890         random orig_rm(rm);
 891         const uint N = 1000000000;
 892         printf("Running cycle test, %u tries:\n", N);
 893         for (uint i = 0; i < N; i++)
 894         {
 895             cur_rm.urand32();
 896             if (memcmp(&cur_rm, &orig_rm, sizeof(cur_rm)) == 0)
 897             {
 898                 printf("Cycle detected!\n");
 899                 return false;
 900             }
 901             if ((i & 0xFFFFFF) == 0xFFFFFF)
 902                 printf("%3.2f%%\n", (double)i / N * 100.0f);
 903         }
 904         printf("OK\n");
 905
 906         return true;
 907     }
 908
 909 } // namespace vogl