Ok, so the name "arctan" was a really poor choice for the first function, since it computes arctan(1/x) rather than arctan(x), but anyway... :-)
/ Marcus Comstedt (ACROSS) (Hail Ilpalazzo!)
Previous text:
2003-02-25 14:59: Subject: Re: Do we have a floatingpoint bug?
// Calculate arctan(1/x) of arbitrary precision
Gmp.mpf arctan(int x, int prec) { // The value of the rest term is 1 / (k * x^k) // so we need to iterate at least until x^k > 2^prec / k // 1 is an upper bound for 1/k, so k > prec/log2(x) is good enough. int kbound = 1+prec/(int)(log((float)x)/log(2.0));
Gmp.mpf res = Gmp.mpf(0,prec); int sign=1, xpow=x; for(int k=1; k<=kbound; k+=2, xpow *= x*x, sign=-sign) res += Gmp.mpq(sign, k * xpow);
return res; };
// Calculate PI of arbitrary precision
Gmp.mpf pi(int prec) { return arctan(5, prec)*16-arctan(239, prec)*4; }
/ Marcus Comstedt (ACROSS) (Hail Ilpalazzo!)