Using Pi as a (very slow) random number generator.

gistfile1.cpp

/*  
    This program implements the BBP algorithm to generate a few hexadecimal
    digits beginning immediately after a given position id, or in other words
    beginning at position id + 1.  On most systems using IEEE 64-bit floating-
    point arithmetic, this code works correctly so long as d is less than
    approximately 1.18 x 10^7.  If 80-bit arithmetic can be employed, this limit
    is significantly higher.  Whatever arithmetic is used, results for a given
    position id can be checked by repeating with id-1 or id+1, and verifying 
    that the hex digits perfectly overlap with an offset of one, except possibly
    for a few trailing digits.  The resulting fractions are typically accurate 
    to at least 11 decimal digits, and to at least 9 hex digits.  
*/

/*  David H. Bailey     2006-09-08 */

/* 
 * Simon Knapp 2013-07-20
 *
 * Was interested in assessing a claim that using the digits of pie would make a
 * good random number generator.
 *
 * downloaded source from http://www.experimentalmath.info/bbp-codes/piqpr8.c
 * and modified to what is here.
 *
 * Was going to test using dieharder, but was not prepared to wait to generate
 * an appropriate sample size... it was suggested somewhere that one would need
 * around 25,000,000 samples to run all tests, and on a ubuntu vm running in
 * in virtual box on top of windows on a dell XPS 17 (i7) it took around 5s to
 * generate 1000 random unsigned ints!
 */

#include <math.h>
#include <iostream>
#include <fstream>
#include <cstdlib>

static double get_fraction(int id) {
  double pid, s1, s2, s3, s4;
  double series (int m, int n);
  void ihex (double x, int m, char c[]);
#define NHX 16

/*  id is the digit position.  Digits generated follow immediately after id. */

  s1 = series (1, id);
  s2 = series (4, id);
  s3 = series (5, id);
  s4 = series (6, id);
  pid = 4. * s1 - 2. * s2 - s3 - s4;
  pid = pid - (int) pid + 1.;
  return pid;
}

double series (int m, int id)

/*  This routine evaluates the series  sum_k 16^(id-k)/(8*k+m) 
    using the modular exponentiation technique. */

{
  int k;
  double ak, p, s, t;
  double expm (double x, double y);
#define eps 1e-17

  s = 0.;

/*  Sum the series up to id. */

  for (k = 0; k < id; k++){
    ak = 8 * k + m;
    p = id - k;
    t = expm (p, ak);
    s = s + t / ak;
    s = s - (int) s;
  }

/*  Compute a few terms where k >= id. */

  for (k = id; k <= id + 100; k++){
    ak = 8 * k + m;
    t = pow (16., (double) (id - k)) / ak;
    if (t < eps) break;
    s = s + t;
    s = s - (int) s;
  }
  return s;
}

double expm (double p, double ak)

/*  expm = 16^p mod ak.  This routine uses the left-to-right binary 
    exponentiation scheme. */

{
  int i, j;
  double p1, pt, r;
#define ntp 25
  static double tp[ntp];
  static int tp1 = 0;

/*  If this is the first call to expm, fill the power of two table tp. */

  if (tp1 == 0) {
    tp1 = 1;
    tp[0] = 1.;

    for (i = 1; i < ntp; i++) tp[i] = 2. * tp[i-1];
  }

  if (ak == 1.) return 0.;

/*  Find the greatest power of two less than or equal to p. */

  for (i = 0; i < ntp; i++) if (tp[i] > p) break;

  pt = tp[i-1];
  p1 = p;
  r = 1.;

/*  Perform binary exponentiation algorithm modulo ak. */

  for (j = 1; j <= i; j++){
    if (p1 >= pt){
      r = 16. * r;
      r = r - (int) (r / ak) * ak;
      p1 = p1 - pt;
    }
    pt = 0.5 * pt;
    if (pt >= 1.){
      r = r * r;
      r = r - (int) (r / ak) * ak;
    }
  }

  return r;
}


/*
 * compile with something like:
 * gcc -O3 -o pirand pi_test.cpp -lm -lstdc++
 *
 * and run like:
 * ./pirand <number of samples> <output file name>
 *
 * e.g.
 * ./pirand 1000 out.txt
 */

int main(int argc, char** argv) {
  int n(atoi(argv[1]));
    unsigned int r;
    double d;
    std::ofstream out(argv[2], std::ios::out);
    // header required for dieharder.
    out << "type: d" << '\n'
        << "count: " << n << '\n'
        << "numbit: 32" << '\n';
    
    for(int i=0; i<n; ++i) {
        d = fabs(get_fraction(i*8));
        r = static_cast<unsigned int>(4294967296.0 * (d - floor(d)));
        out << r << '\n';
    }
    out.close();
    return 0;
}