piqpr8.c - Generating Hex Digits of PI

A C source code is provided to generate hexadecimal digits of the fractional part of PI (0.1415927...) using the BBP (Bailey–Borwein–Plouffe) formula.

If you are interested in generating hexadecimal digits of the fractional part of PI (0.1415927...) with a C program, here is the source code, piqpr8.c, provided by David H. Bailey using the BBP (Bailey–Borwein–Plouffe) formula:

/*  
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 */

#include <stdio.h>
#include <math.h>

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

/*  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.;
  ihex (pid, NHX, chx);
printf(" position = %i\n fraction = %.15f \n hex digits = %10.10s\n",
  id, pid, chx);
}

void ihex (double x, int nhx, char chx[])

/*  This returns, in chx, the first nhx hex digits of the fraction 
of x. */

{
  int i;
  double y;
  char hx[] = "0123456789ABCDEF";

  y = fabs (x);

  for (i = 0; i < nhx; i++){
    y = 16. * (y - floor (y));
    chx[i] = hx[(int) y];
  }
}

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, eps, 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;
}

For more details on piqpr8.c, please see BBP Code Directory.

Of course, you can also read the Wikipedia Bailey–Borwein–Plouffe formula page.

Last update: 2015.

Table of Contents

 About This Book

Blowfish Cipher Algorithm

 Blowfish Cipher Encryption Algorithm

 Key Schedule (Sub-Keys and S-Boxes) Algorithm

 Efficient Form of the Blowfish Algorithm

 Blowfish Cipher Decryption Algorithm

 Proof of Blowfish Cipher Algorithm

 Blowfish Cipher Test Vectors

 First 8336 Fractional Hex Digits of PI

piqpr8.c - Generating Hex Digits of PI

 Perl Crypt::Blowfish Module

 Perl Crypt::ECB Perl Module

 Perl Crypt::CBC Module

 Perl Crypt::CFB Perl Module

 OpenSSL "enc -bf-ecb" for Blowfish/ECB Encryption

 OpenSSL "enc -bf-cbc" for Blowfish/CBC Encryption

 OpenSSL "enc -bf-cfb" for Blowfish/CFB Encryption

 OpenSSL "enc -bf-ofb" for Blowfish/OFB Encryption

 References

 PDF Printing Version