Question

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.

Was it helpful?

Solution

This line here looks like a likely culprit:

char* buffer = (char*) malloc(sizeof(char));

Here you allocate space for only a single character.

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top