Numerical problems with the GNU MP Bignum library

Question

I make a series of crunches that require high precision arithmetic. I use the GNU MP library and according to the GMP manual :

"A floating point or float for short is an arbitrary mantissa precision with a limited precision.

Although the mantissa should have arbitrary precision, I still run into accuracy problems. Instead of porting you with my actual code, here is an approximate minimum working example that illustrates my problem. The code computes 9.3 ^ 15, 9.8 ^ 15 and (9.3 * 9.8) ^ 15. On my machine, the values ​​(9.3 ^ 15) * (9.8 ^ 15) and (9.3 * 9.8) ^ 15 begin to differ from the 16th digit and in in this case lead to an error (around) 4.94 * 10 ^ 13.

Any help would be greatly appreciated. The code is 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 see is due to the fact that 9.3 * 9.8 is calculated approximately . Please change the literals as instances of mpf_class:

 mpf_class a, b; a = 9.3; b = 9.8; // ... x *= a; y *= b; z *= a * b; 

If you want infinite precision, consider 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 

Answer 2

The problem is that you are not explicitly setting the accuracy, so you get the standard accuracy, which is usually (I think) 64 bits, so the results are different in the last bit (s) due to different rounding in different calculation methods. This is approximately 20 digits of the general prefix (the difference can increase significantly with large calculations). If you set higher accuracy,

 #include <gmp.h> #include <gmpxx.h> #include <iostream> #include <iomanip> int main() { mpf_class x(1.0,200), y(1.0,200), z(1.0,200), a("9.3",200), b("9.8",200), c(0,200); c = a*b; for (int i = 0; i < 15; i++) { x *= a; y *= b; z *= c; } std::cout << z << "\n" << (x*y) << std::endl; std::cout << z - x*y << std::endl; return 0; } 

200 bits here, you will get a more accurate result:

 $ ./a.out 2.48677e+29 2.48677e+29 -4.80637e-49 

therefore, the general prefix is ​​about 80 decimal digits or almost 256 bits (the smallest multiple of 64 is greater than 199).

With an accuracy of 2000, the difference is -2.78942e-588 with line constructors, 0 if they are initialized from double (but then, of course, the initial precision is limited to 53 bits, so that just means that both paths accumulate the error in the same way).