GNU MP low precision while using mpf_pow function

Question

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?

I think that you're missing a '1' at the beginning of your result (otherwise it disagrees with WolframAlpha by **a lot** more than "the 15th digit"). — , Jan 11 '22 at 16:47
12.3 is not exactly representable in binary floating point. The closest representable IEEE-754 binary64 approximation to it (which is likely what your C implementation uses) is good to about 15 decimal digits. I speculate, then, that Wolfram is being more careful than you are about representing that value (GMP should be able to do it exactly). You get a different result, then, because you are performing a different computation. — John Bollinger, Jan 11 '22 at 16:51
Yes! You are correct @bbbbbbbbb, thanks for the heads up. It was correct in the original [answer](https://stackoverflow.com/questions/70520908/power-function-for-mpf-class-numbers-in-gmpxx/70669988#70669988) though... yet another copy & paste mistake. — Arc, Jan 11 '22 at 16:51
BTW, I believe that 12.3 ^ 123 doesn't require `mpf_set_default_prec` of more than 200. — , Jan 11 '22 at 16:52
@bbbbbbbbb, yes, I guess it doesn't, but I went for the 'exaggeration method' on that one, of course it wouldn't work. — Arc, Jan 11 '22 at 16:54
@bbbbbbbbb: We can see the result, around 1.14•10^134, has a bit in the 2^445 position, so it requires 446 bits just to represent the integer portion, many more if you want the fraction part precisely too. — Eric Postpischil, Jan 11 '22 at 17:01
@JohnBollinger, Mathematica actually uses GMP internally, see [here](https://www.wolfram.com/legal/third-party-licenses/gmp.html), so which one is correct? — Arc, Jan 11 '22 at 17:01
@Arc: Mathematica does not pass “12.3” through `double` first. — Eric Postpischil, Jan 11 '22 at 17:02
@EricPostpischil, yes! That seems to the correct answer! That makes sense, since it does not use *all* of GMP functions, but has its own implementations for many functions. — Arc, Jan 11 '22 at 17:05

John Bollinger · Accepted Answer · 2022-01-11T19:41:48.077

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.

Excellent! A lesson learned: using the GMP string assigment is far more accurate than the assigment with `double`, in fact it must convert the string to mpf all the way up to the precision specified. I'll update the original [answer](https://stackoverflow.com/questions/70520908/power-function-for-mpf-class-numbers-in-gmpxx) with due reference to this one. Thanks! — Arc, Jan 11 '22 at 18:31
GMP cannot represent 12.3 exactly, but it can represent a much closer approximation than double. — Marc Glisse, Jan 11 '22 at 18:41
@MarcGlisse, what I mean is that I assume the way GMP works is to perform the conversion from decimal string to binary float `mpf_t` up to the precision previously specified for that `mpf_t`, but please correct me if I'm wrong. — Arc, Jan 11 '22 at 19:00
That's true, but that's very different from saying "GMP can represent 12.3 (decimal) exactly" and "you get an exact representation". There exist libraries that use a decimal representation, just not GMP. — Marc Glisse, Jan 11 '22 at 19:03

score 0 · Answer 2 · answered Jan 11 '22 at 16:50

0

This gives a result similar to WolframAlpha's:

from decimal import Decimal
from decimal import getcontext

getcontext().prec = 200

print(Decimal('12.3') ** 123)

So you must be doing something wrong in your GMP configuration.

answered Jan 11 '22 at 16:50

This does not seem to attempt to answer the question. – John Bollinger Jan 11 '22 at 16:52
@JohnBollinger: First off, if you want to be 100% strict, then it does in fact answer the question `is WolframAlpha giving an incorrect result?`. Second, I don't have any way to properly post it in a comment. – Jan 11 '22 at 16:53
What language is it, Mathematica? I only have access to WolframAlpha. And to what result does it agree? Can you post the actual result? – Arc Jan 11 '22 at 16:55
@Arc: Python. `114374367934617190099880295228066276746218078451850229775887975052369504785666896446606568365201542169649974727730628842345343196581134.89591994282087444983721209947664895835902379607854904194900780722`. – Jan 11 '22 at 16:59
So Python agrees to the last digit you posted with WolframAlpha. – Arc Jan 11 '22 at 17:04
@Arc: Well, I didn't go too deep into comparing those 200 digits, but the first 15 definitely looked the same. Though I suppose that one could argue that perhaps WolframAlpha uses Python, and both of them are wrong. – Jan 11 '22 at 17:05
I just paste them into a text editor, they match. But I guess that @EricPostpischil already solved the issue. – Arc Jan 11 '22 at 17:07
Let us [continue this discussion in chat](https://chat.stackoverflow.com/rooms/240965/discussion-between-arc-and-bbbbbbbbb), if anything else has to be said. I've already been [instructed](https://meta.stackoverflow.com/questions/415189/very-important-discussion-deleted-from-comments-a-possible-bug-in-glibcs-mallo) by moderators not to make long conversations on comments. – Arc Jan 11 '22 at 17:18

GNU MP low precision while using mpf_pow function

2 Answers2

Linked