Потеря точности (с использованием GMP) при расчете пи с последовательным bbp

Question

Я пытаюсь реализовать этот алгоритм в C (вместо Python) с использованием GMP. Их код на Python довольно короткий; тем не менее, реализация C намного длиннее, потому что мы не можем выполнять кастинг опор GMPY.

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

Выход становится неверным после 3.243F6A8885A308D313198A.... Что происходит?

РЕДАКТИРОВАТЬ: я изменил метод f_mod(), чтобы он вычитал два mpqs вместо mpfs. Это немного увеличило точность, но все же недостаточно.

РЕДАКТИРОВАТЬ 2: Я попытался запустить эту программу на моем другом ПК, и он волшебным образом работал! Не уверен, что отличается от того факта, что ноутбук, на котором я это написал, был 32-разрядным, а рабочий стол, на котором он работал, был 64-разрядным. Это также может быть ошибкой в ​​реализации GMP (версия, упакованная в Debian Wheezy). Однако я подозреваю, что с этой последовательной формулой BBP что-то не так. Я думаю, что я попытаюсь восстановить это.

Answer 1

Эта строка выглядит как вероятный виновник:

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

Здесь вы выделяете место только для одного символа.