MPFR Accuracy and precision problem, binding to pyton with ctypes

I am trying to call the "fpow" C function from python with "fpow" wrapper. It works but unfortunately both precision and accuracy is lacking. I guess that my use of mpfr is incorrect. How can I solve this?

numeric.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gmp.h>
#include <mpfr.h>
#include "numeric.h"

char *fpow(char *x, char *y, unsigned long long prec)
{
    mpfr_t a, b;
    mpfr_init2(a, prec);
    mpfr_init2(b, prec);

    mpfr_set_str(a, x, 10, 0);
    mpfr_set_str(b, y, 10, 0);

    mpfr_pow(a, a, b, 1);

    long dot;
    char *str = malloc(sizeof(char) * (prec + 1));
    mpfr_get_str(str, &dot, 10, prec, a, 1);

    for(long j = dot; j < dot; j--)
        str[j + 1] = str[j];

    str[dot] = '.';

    mpfr_clear(a);
    mpfr_clear(b);
    mpfr_free_cache();

    return str;
}

basics.py

import ctypes as c


cfunc = c.CDLL('numeric.so')


def fpow(base: str, exponent: str, precision: int) -> str:
    '''base to the exponent with the specified precision'''

    cfunc.fpow.argtypes = (c.c_char_p, c.c_char_p, c.c_ulonglong)
    cfunc.fpow.restype = c.c_char_p

    return (
        cfunc.fpow(
            c.create_string_buffer(base.encode('utf-8')),
            c.create_string_buffer(exponent.encode('utf-8')),
            precision
        )).decode('utf-8')

test.py

import Numer.basics as nm

print(nm.fpow('5.1', '3', 100))

Output of the test.py is "132.509999999999999999999999997382748843093935872477183435247383158639422617852687835693359375000000"

Solution

Your for loop has a bug. It starts with j equal to dot, so it will execute zero times. And, it is decrementing j.

Strings in C need an extra 0x00 char at the end to denote "end of string".

The malloc allocates space for prec + 1 extra char. This is enough for the mpfr_pow output [with 0x00 EOS string terminator].

But this does not allocate space for the . we are trying to add. So, we need to increase the amount allocated by malloc.

The [redacted] raw output of mpfr_pow is:

132650999999999999999999999999738274884309393587247718343524738315 ...

But, the output of your function is:

132.50999999999999999999999999738274884309393587247718343524738315 ...

Notice that you're trashing the [first] 6 with .

Here's a refactored version of the C code. I added main to be able to test independently of the python code:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gmp.h>
#include <mpfr.h>
//#include "numeric.h"

char *
fpow(char *x, char *y, unsigned long long prec)
{
    mpfr_t a, b;

    mpfr_init2(a, prec);
    mpfr_init2(b, prec);

    mpfr_set_str(a, x, 10, 0);
    mpfr_set_str(b, y, 10, 0);

    mpfr_pow(a, a, b, 1);

    long dot;
// NOTE/BUG: allocates enough space for 0x00 EOS char but _not_ the added
// '.' character
#if 0
    char *str = malloc(sizeof(char) * (prec + 1));
#else
    char *str = malloc(sizeof(char) * (prec + 1 + 1));
#endif

    mpfr_get_str(str, &dot, 10, prec, a, 1);

// original code ...
#if 0
    for (long j = dot; j < dot; j--)
        str[j + 1] = str[j];
    str[dot] = '.';
// refactored code ...
#else
    int prev = str[dot];
    int next;
    char *cur;
    for (cur = &str[dot + 1];  *cur != 0;  ++cur) {
        next = *cur;
        *cur = prev;
        prev = next;
    }
    *cur++ = prev;
    *cur = 0;
    str[dot] = '.';
#endif

    mpfr_clear(a);
    mpfr_clear(b);
    mpfr_free_cache();

    return str;
}

int
main(void)
{

    char *ret = fpow("5.1","3",100);
    printf("%s\n",ret);

    return 0;
}

UPDATE:

I have 2 questions. First is there a reason to not fix the for loop bug as such: for (long j = dot; j < prec; j--) str[j + 1] = str[j]; str[dot] = '.'; This seems to get the work done aswell.

It didn't work for me. Of course, there are different ways to do the insert. Mine [which works] is but one. Feel free to experiment. You could isolate the insertion code into a separate function [and try different functions].

The second question is that as you can probably try and see "prec" doesn't really work e.g. when i put prec=20 it doesn't give 20 digits accurately. – aras edeş

That's because for mpfr_init2, the second argument is not the number of decimal digits. It is the number of bits.

So, roughly, to get a given number of digits of precision, we need approximately 4x the number of bits (e.g. to represent the digit 9 we need 4 bits).

So, change:

mpfr_init2(a, prec);
mpfr_init2(b, prec);

Into:

mpfr_init2(a, prec * 4);
mpfr_init2(b, prec * 4);

Note that we could calculate the needed precision exactly but I think that's overkill in the given use case.

With that change, the redacted output is:

132.65099999999999999999999999999999999999999999999999999999999999 ...