Numerical Problems with GNU MP Bignum Library

Question

I'm doing some number crunching, which needs high precision arithmetic. I'm using the GNU MP library, and according to the GMP manual:

"Floating point number or Float for short, is an arbitrary precision mantissa with a limited precision exponent."

Although the mantissa is supposed to have an arbitrary precision, I'm still running into precision problems. Rather than bore you with my actual code, here is a near-minimal working example that illustrates my problem. The code calculates 9.3^15, 9.8^15 and (9.3*9.8)^15. On my machine, the values of (9.3^15)*(9.8^15) and (9.3*9.8)^15 start to differ from the 16-th digit onwards, in this case leading to an error of (around) 4.94*10^13.

Any help will be greatly appreciated. Code below.

#include <gmp.h>
#include <gmpxx.h>

#include <iostream>
#include <iomanip>

int main()
{
    mpf_class x, y, z;
    x = y = z = 1.0;

    for (int i = 0; i < 15; i++)
    {
        x *= 9.3;
        y *= 9.8;
        z *= 9.3*9.8;
    }

    std::cout << z - x*y << std::endl;

    return 0;
}

Answer 1

The problem you are seeing is due to the fact that 9.3 * 9.8 is calculated approximately. Please change the literals to be instances of the mpf_class:

mpf_class a, b;
a = 9.3;
b = 9.8;

// ...

x *= a;
y *= b;
z *= a * b;

If you need infinite precision, consider using rational numbers instead:

#include <gmp.h>
#include <gmpxx.h>

#include <iostream>
#include <iomanip>

int main()
{
    mpq_class x(1), y(1), z(1), a(93, 10), b(98, 10);

    for (int i = 0; i < 15; i++)
    {
        x *= a;
        y *= b;
        z *= (a * b);
    }

    std::cout << z - x*y << std::endl << z << std::endl;

    return 0;
}

prints

0
7589015305950762920038660273144124106674963183136666693/30517578125000000000000000