I'm trying to implement this algorithm in C (instead of Python) using GMP. Their Python code is quite short; however, the C implementation is a lot longer because we can't do the casting GMPY supports.
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
void f_mod(mpq_t input) {
mpf_t numerator;
mpf_init(numerator);
mpf_set_z(numerator, mpq_numref(input));
mpf_t denominator;
mpf_init(denominator);
mpf_set_z(denominator, mpq_denref(input));
mpf_t div;
mpf_init(div);
mpf_t integer_part;
mpf_init(integer_part);
// Integer component
mpq_t floored;
mpq_init(floored);
mpf_div(div, numerator, denominator);
mpf_floor(integer_part, div);
mpq_set_f(floored, integer_part);
// Get fractional part
mpq_sub(input, input, floored);
}
int main() {
unsigned long end_val = 1000;
mpq_t x;
mpq_init(x);
mpq_set_ui(x, 0, 1);
mpq_t p;
mpq_init(p);
mpq_t res;
mpq_init(res);
mpq_t sixteen;
mpq_init(sixteen);
mpq_set_ui(sixteen, 16, 1);
unsigned long n = 1;
printf("\n");
short* results = (short*) malloc(end_val * sizeof(short));
while (n < end_val)
{
mpq_set_ui(p, (120*n-89)*n+16, (((512*n-1024)*n+712)*n-206)*n+21);
mpq_mul(res, x, sixteen);
mpq_add(res, res, p);
f_mod(res);
mpq_set(x, res);
mpq_mul(res, res, sixteen);
results[n] = (short) mpq_get_d(res);
n++;
}
char* buffer = (char*) malloc(sizeof(char));
for (unsigned long a = 0; a < end_val; a++)
{
snprintf(buffer, 16, "%x", results[a]);
printf("%s", buffer);
}
return 0;
}
The output becomes incorrect after 3.243F6A8885A308D313198A.... What's happening?
EDIT: I changed the f_mod() method so it subtracts two mpqs instead of mpfs. This increased accuracy a little bit, but still not enough.
EDIT 2: I tried running this program on my other PC, and it magically worked! Not sure what's different other than the fact the laptop I wrote this on was 32-bit and the desktop on which it worked was 64-bit. It also might be a bug in the GMP implementation (version packaged into Debian Wheezy). However, I do suspect there is something wrong with that sequential BBP formula. I think I'll try to re-derive it.
This line here looks like a likely culprit:
char* buffer = (char*) malloc(sizeof(char));
Here you allocate space for only a single character.