X-Git-Url: https://ginac.de/CLN/cln.git//cln.git?a=blobdiff_plain;f=src%2Ffloat%2Ftranscendental%2Fcl_LF_pi.cc;h=3b537ad52e5929fbdddc90af65be9e3bbfc7cd7f;hb=3af2cde18b3aabed4c808b0113daa81c2263b0bd;hp=6660f62911418e093f9d147cb7c7a453d9bed22a;hpb=38f5f03b3f6d25840a83294a08b136a44f79de4f;p=cln.git

diff --git a/src/float/transcendental/cl_LF_pi.cc b/src/float/transcendental/cl_LF_pi.cc
index 6660f62..3b537ad 100644
--- a/src/float/transcendental/cl_LF_pi.cc
+++ b/src/float/transcendental/cl_LF_pi.cc
@@ -1,26 +1,26 @@
 // pi().
 
 // General includes.
-#include "cl_sysdep.h"
+#include "base/cl_sysdep.h"
 
 // Specification.
-#include "cl_F_tran.h"
+#include "float/transcendental/cl_F_tran.h"
 
 
 // Implementation.
 
 #include "cln/lfloat.h"
-#include "cl_LF_tran.h"
-#include "cl_LF.h"
+#include "float/transcendental/cl_LF_tran.h"
+#include "float/lfloat/cl_LF.h"
 #include "cln/integer.h"
-#include "cl_alloca.h"
+#include "base/cl_alloca.h"
 
 namespace cln {
 
 ALL_cl_LF_OPERATIONS_SAME_PRECISION()
 
 // For the next algorithms, I warmly recommend
-// [Jörg Arndt: hfloat documentation, august 1996,
+// [Jörg Arndt: hfloat documentation, august 1996,
 //  http://www.tu-chemnitz.de/~arndt/ hfdoc.dvi
 //  But beware of the typos in his formulas! ]
 
@@ -31,7 +31,7 @@ const cl_LF compute_pi_brent_salamin (uintC len)
 	//  functions. J. ACM 23(1976), 242-251.]
 	// [Jonathan M. Borwein, Peter B. Borwein: Pi and the AGM.
 	//  Wiley 1987. Algorithm 2.2, p. 48.]
-	// [Jörg Arndt, formula (4.51)-(4.52).]
+	// [Jörg Arndt, formula (4.51)-(4.52).]
 	//    pi = AGM(1,1/sqrt(2))^2 * 2/(1 - sum(k=0..infty, 2^k c_k^2)).
 	// where the AGM iteration reads
 	//    a_0 := 1, b_0 := 1/sqrt(2).
@@ -61,9 +61,9 @@ const cl_LF compute_pi_brent_salamin (uintC len)
 	//   )
 	//   (/ (expt a 2) t)
 	// )
-	var uintC actuallen = len + 1; // 1 Schutz-Digit
+	var uintC actuallen = len + 1; // 1 guard digit
 	var uintE uexp_limit = LF_exp_mid - intDsize*len;
-	// Ein Long-Float ist genau dann betragsmäßig <2^-n, wenn
+	// Ein Long-Float ist genau dann betragsmäßig <2^-n, wenn
 	// sein Exponent < LF_exp_mid-n = uexp_limit ist.
 	var cl_LF a = cl_I_to_LF(1,actuallen);
 	var cl_LF b = sqrt(scale_float(a,-1));
@@ -79,14 +79,14 @@ const cl_LF compute_pi_brent_salamin (uintC len)
 		k++;
 	}
 	var cl_LF pires = square(a)/t; // a^2/t
-	return shorten(pires,len); // verkürzen und fertig
+	return shorten(pires,len); // verkürzen und fertig
 }
 // Bit complexity (N := len): O(log(N)*M(N)).
 
 const cl_LF compute_pi_brent_salamin_quartic (uintC len)
 {
 	// See [Borwein, Borwein, section 1.4, exercise 3, p. 17].
-	// See also [Jörg Arndt], formulas (4.52) and (3.30)[wrong!].
+	// See also [Jörg Arndt], formulas (4.52) and (3.30)[wrong!].
 	// As above, we are using the formula
 	//    pi = AGM(1,1/sqrt(2))^2 * 2/(1 - sum(k=0..infty, 2^k c_k^2)).
 	// where the AGM iteration reads
@@ -111,7 +111,7 @@ const cl_LF compute_pi_brent_salamin_quartic (uintC len)
 	//   = 2^(k+1) * [wa_k^4 - ((wa_k^2+wb_k^2)/2)^2].
 	// Hence,
 	//   pi = AGM(1,1/sqrt(2))^2 * 1/(1/2 - sum(k even, 2^k*[....])).
-	var uintC actuallen = len + 1; // 1 Schutz-Digit
+	var uintC actuallen = len + 1; // 1 guard digit
 	var uintE uexp_limit = LF_exp_mid - intDsize*len;
 	var cl_LF one = cl_I_to_LF(1,actuallen);
 	var cl_LF a = one;
@@ -135,7 +135,7 @@ const cl_LF compute_pi_brent_salamin_quartic (uintC len)
 		k += 2;
 	}
 	var cl_LF pires = square(a)/t;
-	return shorten(pires,len); // verkürzen und fertig
+	return shorten(pires,len); // verkürzen und fertig
 }
 // Bit complexity (N := len): O(log(N)*M(N)).
 
@@ -146,7 +146,7 @@ const cl_LF compute_pi_ramanujan_163 (uintC len)
 	// mit J = -53360^3 = - (2^4 5 23 29)^3
 	//     A = 163096908 = 2^2 3 13 1045493
 	//     B = 6541681608 = 2^3 3^3 7 11 19 127 163
-	// See [Jörg Arndt], formulas (4.27)-(4.30).
+	// See [Jörg Arndt], formulas (4.27)-(4.30).
 	// This is also the formula used in Pari.
 	// The absolute value of the n-th summand is approximately
 	//   |J|^-n * n^(-1/2) * B/(2*pi^(3/2)),
@@ -154,7 +154,7 @@ const cl_LF compute_pi_ramanujan_163 (uintC len)
 	// in precision.
 	// The sum is best evaluated using fixed-point arithmetic,
 	// so that the precision is reduced for the later summands.
-	var uintC actuallen = len + 4; // 4 Schutz-Digits
+	var uintC actuallen = len + 4; // 4 guard digits
 	var sintC scale = intDsize*actuallen;
 	static const cl_I A = "163096908";
 	static const cl_I B = "6541681608";
@@ -177,16 +177,10 @@ const cl_LF compute_pi_ramanujan_163 (uintC len)
 	var cl_LF fsum = scale_float(cl_I_to_LF(sum,actuallen),-scale);
 	static const cl_I J3 = "262537412640768000"; // -1728*J
 	var cl_LF pires = sqrt(cl_I_to_LF(J3,actuallen)) / fsum;
-	return shorten(pires,len); // verkürzen und fertig
+	return shorten(pires,len); // verkürzen und fertig
 }
 // Bit complexity (N := len): O(N^2).
 
-#if defined(__mips__) && !defined(__GNUC__) // workaround SGI CC bug
-#define A A_fast
-#define B B_fast
-#define J3 J3_fast
-#endif
-
 const cl_LF compute_pi_ramanujan_163_fast (uintC len)
 {
 	// Same formula as above, using a binary splitting evaluation.
@@ -216,7 +210,7 @@ const cl_LF compute_pi_ramanujan_163_fast (uintC len)
 			: cl_pqa_series_stream (rational_series_stream::computenext),
 			  n (0) {}
 	} series;
-	var uintC actuallen = len + 4; // 4 Schutz-Digits
+	var uintC actuallen = len + 4; // 4 guard digits
 	static const cl_I A = "163096908";
 	static const cl_I B = "6541681608";
 	static const cl_I J1 = "10939058860032000"; // 72*abs(J)
@@ -230,10 +224,10 @@ const cl_LF compute_pi_ramanujan_163_fast (uintC len)
 	var uintC N = (n_slope*actuallen)/32 + 1;
 	// N > intDsize*log(2)/log(|J|) * actuallen, hence
 	// |J|^-N < 2^(-intDsize*actuallen).
-	var cl_LF fsum = eval_rational_series(N,series,actuallen);
+	var cl_LF fsum = eval_rational_series<false>(N,series,actuallen,actuallen);
 	static const cl_I J3 = "262537412640768000"; // -1728*J
 	var cl_LF pires = sqrt(cl_I_to_LF(J3,actuallen)) / fsum;
-	return shorten(pires,len); // verkürzen und fertig
+	return shorten(pires,len); // verkürzen und fertig
 }
 // Bit complexity (N := len): O(log(N)^2*M(N)).
 
@@ -260,22 +254,22 @@ const cl_LF compute_pi_ramanujan_163_fast (uintC len)
 
 const cl_LF pi (uintC len)
 {
-	var uintC oldlen = TheLfloat(cl_LF_pi)->len; // vorhandene Länge
+	var uintC oldlen = TheLfloat(cl_LF_pi())->len; // vorhandene Länge
 	if (len < oldlen)
-		return shorten(cl_LF_pi,len);
+		return shorten(cl_LF_pi(),len);
 	if (len == oldlen)
-		return cl_LF_pi;
+		return cl_LF_pi();
 
-	// TheLfloat(cl_LF_pi)->len um mindestens einen konstanten Faktor
-	// > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird:
+	// TheLfloat(cl_LF_pi())->len um mindestens einen konstanten Faktor
+	// > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird:
 	var uintC newlen = len;
 	oldlen += floor(oldlen,2); // oldlen * 3/2
 	if (newlen < oldlen)
 		newlen = oldlen;
 
-	// gewünschte > vorhandene Länge -> muß nachberechnen:
-	cl_LF_pi = compute_pi_ramanujan_163_fast(newlen);
-	return (len < newlen ? shorten(cl_LF_pi,len) : cl_LF_pi);
+	// gewünschte > vorhandene Länge -> muß nachberechnen:
+	cl_LF_pi() = compute_pi_ramanujan_163_fast(newlen);
+	return (len < newlen ? shorten(cl_LF_pi(),len) : cl_LF_pi());
 }
 
 }  // namespace cln