BBP

/* This program employs the recently discovered digit extraction scheme
to produce hex digits of pi. This code is valid up to ic = 2^24 on
systems with IEEE arithmetic. */

/* David H. Bailey 960429 */

#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 ic = 1000000;
    #define NHX 16
    char chx[NHX];

    /* ic is the hex digit position -- output begins at position ic + 1. */

    s1 = series (1, ic);
    s2 = series (4, ic);
    s3 = series (5, ic);
    s4 = series (6, ic);
    pid = 4. * s1 - 2. * s2 - s3 - s4;
    pid = pid - (int) pid + 1.;
    ihex (pid, NHX, chx);
    printf ("Pi hex digit computation\n");
    printf ("position = %i + 1\n %20.15f\n %12.12s\n", ic, 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 ic)

/* This routine evaluates the series sum_k 16^(ic-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 ic. */

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

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

    for (k = ic; k <= ic + 100; k++){
        ak = 8 * k + m;
        t = pow (16., (double) (ic - 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. It is valid for ak <= 2^24. */

{
    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;
}