1

I am trying to implement the BBP formula. I wrote the following quick and dirty test code, based on the paper by Paul H. Bailey (that I found here):

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
#include <limits.h>

long bin_exp_mod_k(int base, int exponent, int div) {
    long expmod = 1;
    long t;

    t = (long) pow(2, floor(log(exponent)/log(2.0)));

    while(42) {
        if(exponent >= t) {
            expmod = (base * expmod) % div;
            exponent -= t;
        }

        t = t/2;

        if(t >= 1) {
            expmod = ((long)pow(expmod, 2)) % div;
        } else
            break;
    }

    return expmod;
}

void bbp_part(mpf_t rop, int base, int index, int digit_index, int num_iter) {
    const int knum = num_iter;
    int k;
    long ldiv;
    mpf_t fbase;
    mpf_t fst_sum;
    mpf_t fst_sum_int;
    mpf_t snd_sum;
    mpf_t powval;
    mpf_t div;

    mpf_init(fst_sum);
    mpf_init(snd_sum);
    mpf_init(fst_sum_int);
    mpf_init(powval);
    mpf_init(div);
    mpf_init(rop);
    mpf_init_set_si(fbase, base);

    for(k = 0;k <= digit_index;k++) {
        ldiv = 8 * k + index;
        mpf_set_si(powval, bin_exp_mod_k(base, digit_index - k, ldiv));
        mpf_set_si(div, ldiv);
        mpf_div(powval, powval, div);
        mpf_add(fst_sum, fst_sum, powval);
    }

    mpf_trunc(fst_sum_int, fst_sum);
    mpf_sub(fst_sum, fst_sum, fst_sum_int);

    for(k = digit_index + 1;k < knum;k++) {
        ldiv = 8 * k + index;
        mpf_set_si(div, ldiv);
        mpf_pow_ui(powval, fbase, digit_index - k);
        mpf_div(powval, powval, div);
        mpf_add(snd_sum, snd_sum, powval);
    }

    mpf_set(rop, fst_sum);
    mpf_add(fst_sum, fst_sum, snd_sum);

    printf("S%i: %.20lf\n", index, mpf_get_d(rop));

    mpf_clear(fst_sum);
    mpf_clear(snd_sum);
    mpf_clear(fbase);
    mpf_clear(fst_sum_int);
    mpf_clear(powval);
    mpf_clear(div);
}

int main(int argc, char **argv) {
    const int base = 16;
    int d;
    int num_iter;
    mpf_t pi_digits;
    mpf_t part1;
    mpf_t part1_int;
    mpf_t part4;
    mpf_t part4_int;
    mpf_t part5;
    mpf_t part5_int;
    mpf_t part6;
    mpf_t part6_int;

    if(argc == 1) {
        return -1;
    }

    d = atoi(argv[1]);

    if(argc == 3) {
        num_iter = atoi(argv[2]);
    } else {
        num_iter = INT_MAX;
    }

    mpf_set_default_prec(128);

    mpf_init(pi_digits);

    mpf_init(part1_int);       
    mpf_init(part4_int);
    mpf_init(part5_int);
    mpf_init(part6_int);

    bbp_part(part1, base, 1, d, num_iter);
    bbp_part(part4, base, 4, d, num_iter);
    bbp_part(part5, base, 5, d, num_iter);
    bbp_part(part6, base, 6, d, num_iter);

    mpf_trunc(part1_int, part1);
    mpf_trunc(part4_int, part4);
    mpf_trunc(part5_int, part5);
    mpf_trunc(part6_int, part6);

    mpf_sub(part1, part1, part1_int);
    mpf_sub(part4, part4, part4_int);
    mpf_sub(part5, part5, part5_int);
    mpf_sub(part6, part6, part6_int);

    mpf_mul_ui(part1, part1, 4L);
    mpf_mul_ui(part4, part4, 2L);

    mpf_set(pi_digits, part1);
    mpf_sub(pi_digits, pi_digits, part4);
    mpf_sub(pi_digits, pi_digits, part5);
    mpf_sub(pi_digits, pi_digits, part6);

    mpf_clear(pi_digits);
    mpf_clear(part1);
    mpf_clear(part4);
    mpf_clear(part5);
    mpf_clear(part6);
    mpf_clear(part1_int);
    mpf_clear(part4_int);
    mpf_clear(part5_int);
    mpf_clear(part6_int);

    return 0;
}

and it seems to be correct, because according to his paper the partial results I get for S1, S4, S5 and S6 are, with a small precision difference, the same ones noted.

However, I can't figure out what I did miss with the "combining the results" part. No matter what I do, I get 0.57... instead of the 0.42... written on the paper.

Can someone see what am I missing here?

Mauren
  • 1,955
  • 2
  • 18
  • 28
  • 1
    Check http://stackoverflow.com/questions/7011184/floating-point-comparison , possibly loss of precision due to using floating point types - although i'm not familiar with mpf_t but some googling suggests it's a floating point type. – Walter Delevich Apr 25 '15 at 18:03
  • @WalterDelevich I'm aware of that, however I have no floating point comparisons around. – Mauren Apr 25 '15 at 18:06
  • 1
    It can hardly be the floating point that creates such a big difference. A floating point eror would be in the orde of 10 **-10 at least. – Paul Ogilvie Apr 25 '15 at 18:40
  • I scanned the paper. It gives the mathematics of the formula but does not give an algorithm. Turning a paper into an algorithm is often hard and most likely you do something wrong. Finding where you misintepreted the formula and wrote the wrong code for it, is not the task of StackExchange. – Paul Ogilvie Apr 25 '15 at 18:43
  • @PaulOgilvie I supposed someone could have run into the same problem, since many people could have implemented it. – Mauren Apr 25 '15 at 18:45

0 Answers0