GNU MP low precision while using mpf_pow function

While writing this answer, I used the mpf_pow function to calculate 12.3 ^ 123, and the result is different from the one given by WolframAlpha (which by the way also uses GMP).

I casted the code to pure C to simplify:

#include <stdio.h>
#include <gmp.h>

int main (void) {
    mpf_t a, c;
    unsigned long int b = 123UL;
    mpf_set_default_prec(100000);
    mpf_inits(a, c, NULL);
    mpf_set_d(a, 12.3);
    mpf_pow_ui(c, a, b);
    gmp_printf("c = %.50Ff\n", c);
    return 0;
}

Which results in

114374367934618002778643226182707594198913258409535335775583252201365538178632825702225459029661601216944929436371688246107986574246790.32099077871758646985223686110515186972735931183764

While WolframAlpha returns

1.14374367934617190099880295228066276746218078451850229775887975052369504785666896446606568365201542169649974727730628842345343196581134895919942820874449837212099476648958359023796078549041949007807220625356526926729664064846685758382803707100766740220839267 × 10^134

which starts to disagree with mpf_pow at the 15th digit.

Am I doing something wrong in the code, is this a limitation of GMP, or is WolframAlpha giving an incorrect result?

Solution

Am I doing something wrong in the code, is this a limitation of GMP, or is WolframAlpha giving an incorrect result?

You are doing something different from what Wolfram is doing (obviously). Your code is not wrong, per se, but it is not doing what you probably think it is doing. Compare the output of this variation:

#include <stdio.h>
#include <gmp.h>

int main (void) {
    mpf_t a, c;
    unsigned long int b = 123UL;
    mpf_set_default_prec(100000);
    mpf_inits(a, c, NULL);
    mpf_set_d(a, 12.3);
    mpf_pow_ui(c, a, b);
    gmp_printf("c  = %.50Ff\n", c);

    putchar('\n');
    mpf_t a1, c1;
    mpf_inits(a1, c1, NULL);
    mpf_set_str(a1, "12.3", 10);
    mpf_pow_ui(c1, a1, b);
    gmp_printf("c' = %.50Ff\n", c1);

    return 0;
}

...

c  = 114374367934618002778643226182707594198913258409535335775583252201365538178632825702225459029661601216944929436371688246107986574246790.32099077871758646985223686110515186972735931183764

c' = 114374367934617190099880295228066276746218078451850229775887975052369504785666896446606568365201542169649974727730628842345343196581134.89591994282087444983721209947664895835902379607855

The difference between the two output values arises because my C implementation and yours represent values of type double in binary floating point, and 12.3 is not exactly representable in binary floating point (see Is floating point math broken?). C provides the closest approximation available, which, assuming 64-bit IEEE 754 representation, matches to about 15 decimal digits of precision. When you initialize a GMP variable with such a value, you get an exact GMP representation of the actual double value, which is only an approximation to 12.3 decimal.

But GMP can represent 12.3 (decimal) to whatever precision you choose.^* You chose a very high precision, so when you use a decimal string to initialize your MP-float variable you get a much closer approximation than when you used a double. Naturally, performing the same operation on those different values produces different results. The GMP result in the latter case appears to agree with the Wolfram result to the full precision in which it is expressed.

Note also that in a general sense, one can also use decimal floating-point, in software or (if you are so equipped) in hardware. The value 12.3 (decimal) can be represented exactly in such a format, but that's not what GMP uses.

^* Or indeed, GMP can represent 12.3 exactly as a MP rational, though that's not what the code above does.