I'm trying to test whether GMP is working with the precision I want or not with the next program:
#include<gmpxx.h>
#include<iostream>
#include<iomanip>
int main()
{
long prec = 1000;
mpf_set_default_prec(prec);
mpf_class a(1);
mpf_class b(mpf_class(1)/sqrt(mpf_class(2)));
mpf_class t(mpf_class(1)/mpf_class(4));
mpf_class p(1);
mpf_class x,y,pi, pi2, sub;
pi2 = mpf_class(3.141592653589793238462643383279502884197169399375105820974944592307816406286);
while ( a - b > mpf_class(1e-250))
{
x = (a + b)/2;
y = sqrt(a*b);
t = t - p*(a-x)*(a-x);
a = x;
b = y;
p *=2;
}
pi = (a+b)*(a+b)/(mpf_class(4)*t);
sub = pi - pi2;
std::cout << std::setprecision(80) << pi << '\n' << pi2 << '\n' << sub << '\n';
return 0;
}
It calculates pi with many digits and then substracts pi2 (pi with 100 digits) from pi (value calculated, it's easy to check that all decimals shown are correct).
Output:
3.141592653589793238462643383279502884197169399375105820974944592307816406286209
3.141592653589793115997963468544185161590576171875
1.2246467991473531772260659322750010582097494459230781640628620899862803482534212e-16
The problem is that it's working with only 16 decimal places for pi2 and I don't know why.
The code was taken from http://ubuntuforums.org/showthread.php?t=459229
I've tried different ways of initiating pi2 but I can't manage to get more than 16 digits.
Use a std::string instead of number literal. Any value of the format aaa.bbb is treated as a double (limited precision) and with f suffix as float (limited precision again). In your case it is probably treated as double.
const std::string number = "3.141592653589793238462643383279502884197169399375105820974944592307816406286";
mpz_class pi2(number);