X-Git-Url: https://ginac.de/CLN/cln.git//cln.git?a=blobdiff_plain;f=src%2Ffloat%2Ftranscendental%2Fcl_LF_pi.cc;h=e82abee7c1eff07fc7331e957cb661418f5ab505;hb=8b3d91dec77438c0fe679b10869ab29e6cdeba58;hp=ef22cb73b8632c2e9bd25f13ef0e9a09d70d8792;hpb=dd9e0f894eec7e2a8cf85078330ddc0a6639090b;p=cln.git

diff --git a/src/float/transcendental/cl_LF_pi.cc b/src/float/transcendental/cl_LF_pi.cc
index ef22cb7..e82abee 100644
--- a/src/float/transcendental/cl_LF_pi.cc
+++ b/src/float/transcendental/cl_LF_pi.cc
@@ -1,4 +1,4 @@
-// cl_pi().
+// pi().
 
 // General includes.
 #include "cl_sysdep.h"
@@ -9,16 +9,18 @@
 
 // Implementation.
 
-#include "cl_lfloat.h"
+#include "cln/lfloat.h"
 #include "cl_LF_tran.h"
 #include "cl_LF.h"
-#include "cl_integer.h"
+#include "cln/integer.h"
 #include "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! ]
 
@@ -29,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).
@@ -60,8 +62,8 @@ const cl_LF compute_pi_brent_salamin (uintC len)
 	//   (/ (expt a 2) t)
 	// )
 	var uintC actuallen = len + 1; // 1 Schutz-Digit
-	var uintL uexp_limit = LF_exp_mid - intDsize*(uintL)len;
-	// Ein Long-Float ist genau dann betragsmäßig <2^-n, wenn
+	var uintE uexp_limit = LF_exp_mid - intDsize*len;
+	// 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));
@@ -76,15 +78,15 @@ const cl_LF compute_pi_brent_salamin (uintC len)
 		a = new_a;
 		k++;
 	}
-	var cl_LF pi = square(a)/t; // a^2/t
-	return shorten(pi,len); // verkürzen und fertig
+	var cl_LF pires = square(a)/t; // a^2/t
+	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
@@ -110,7 +112,7 @@ const cl_LF compute_pi_brent_salamin_quartic (uintC len)
 	// 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 uintL uexp_limit = LF_exp_mid - intDsize*(uintL)len;
+	var uintE uexp_limit = LF_exp_mid - intDsize*len;
 	var cl_LF one = cl_I_to_LF(1,actuallen);
 	var cl_LF a = one;
 	var cl_LF wa = one;
@@ -132,8 +134,8 @@ const cl_LF compute_pi_brent_salamin_quartic (uintC len)
 		b = new_b; wb = new_wb;
 		k += 2;
 	}
-	var cl_LF pi = square(a)/t;
-	return shorten(pi,len); // verkürzen und fertig
+	var cl_LF pires = square(a)/t;
+	return shorten(pires,len); // verkürzen und fertig
 }
 // Bit complexity (N := len): O(log(N)*M(N)).
 
@@ -144,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)),
@@ -152,8 +154,8 @@ 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 uintL actuallen = len + 4; // 4 Schutz-Digits
-	var sintL scale = intDsize*actuallen;
+	var uintC actuallen = len + 4; // 4 Schutz-Digits
+	var sintC scale = intDsize*actuallen;
 	static const cl_I A = "163096908";
 	static const cl_I B = "6541681608";
 	//static const cl_I J1 = "10939058860032000"; // 72*abs(J)
@@ -174,8 +176,8 @@ 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 pi = sqrt(cl_I_to_LF(J3,actuallen)) / fsum;
-	return shorten(pi,len); // verkürzen und fertig
+	var cl_LF pires = sqrt(cl_I_to_LF(J3,actuallen)) / fsum;
+	return shorten(pires,len); // verkürzen und fertig
 }
 // Bit complexity (N := len): O(N^2).
 
@@ -189,7 +191,32 @@ const cl_LF compute_pi_ramanujan_163_fast (uintC len)
 {
 	// Same formula as above, using a binary splitting evaluation.
 	// See [Borwein, Borwein, section 10.2.3].
-	var uintL actuallen = len + 4; // 4 Schutz-Digits
+	struct rational_series_stream : cl_pqa_series_stream {
+		uintC n;
+		static cl_pqa_series_term computenext (cl_pqa_series_stream& thisss)
+		{
+			static const cl_I A = "163096908";
+			static const cl_I B = "6541681608";
+			static const cl_I J1 = "10939058860032000"; // 72*abs(J)
+			var rational_series_stream& thiss = (rational_series_stream&)thisss;
+			var uintC n = thiss.n;
+			var cl_pqa_series_term result;
+			if (n==0) {
+				result.p = 1;
+				result.q = 1;
+			} else {
+				result.p = -((cl_I)(6*n-5)*(cl_I)(2*n-1)*(cl_I)(6*n-1));
+				result.q = (cl_I)n*(cl_I)n*(cl_I)n*J1;
+			}
+			result.a = A+n*B;
+			thiss.n = n+1;
+			return result;
+		}
+		rational_series_stream ()
+			: cl_pqa_series_stream (rational_series_stream::computenext),
+			  n (0) {}
+	} series;
+	var uintC actuallen = len + 4; // 4 Schutz-Digits
 	static const cl_I A = "163096908";
 	static const cl_I B = "6541681608";
 	static const cl_I J1 = "10939058860032000"; // 72*abs(J)
@@ -198,40 +225,15 @@ const cl_LF compute_pi_ramanujan_163_fast (uintC len)
 	//   a(n) = A+n*B, b(n) = 1,
 	//   p(n) = -(6n-5)(2n-1)(6n-1) for n>0,
 	//   q(n) = 72*|J|*n^3 for n>0.
-	var const uintL n_slope = (uintL)(intDsize*32*0.02122673)+1;
+	var const uintC n_slope = (uintC)(intDsize*32*0.02122673)+1;
 	// n_slope >= 32*intDsize*log(2)/log(|J|), normally n_slope = 22.
-	var uintL N = (n_slope*actuallen)/32 + 1;
+	var uintC N = (n_slope*actuallen)/32 + 1;
 	// N > intDsize*log(2)/log(|J|) * actuallen, hence
 	// |J|^-N < 2^(-intDsize*actuallen).
-	CL_ALLOCA_STACK;
-	var cl_I* av = (cl_I*) cl_alloca(N*sizeof(cl_I));
-	var cl_I* pv = (cl_I*) cl_alloca(N*sizeof(cl_I));
-	var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
-	var uintL* qsv = (uintL*) cl_alloca(N*sizeof(uintL));
-	var uintL n;
-	for (n = 0; n < N; n++) {
-		init1(cl_I, av[n]) (A+n*B);
-		if (n==0) {
-			init1(cl_I, pv[n]) (1);
-			init1(cl_I, qv[n]) (1);
-		} else {
-			init1(cl_I, pv[n]) (-((cl_I)(6*n-5)*(cl_I)(2*n-1)*(cl_I)(6*n-1)));
-			init1(cl_I, qv[n]) ((cl_I)n*(cl_I)n*(cl_I)n*J1);
-		}
-	}
-	var cl_pqa_series series;
-	series.av = av;
-	series.pv = pv; series.qv = qv;
-	series.qsv = (len >= 35 ? qsv : 0); // 5% speedup for large len's
-	var cl_LF fsum = eval_rational_series(N,series,actuallen);
-	for (n = 0; n < N; n++) {
-		av[n].~cl_I();
-		pv[n].~cl_I();
-		qv[n].~cl_I();
-	}
+	var cl_LF fsum = eval_rational_series(N,series,actuallen,actuallen);
 	static const cl_I J3 = "262537412640768000"; // -1728*J
-	var cl_LF pi = sqrt(cl_I_to_LF(J3,actuallen)) / fsum;
-	return shorten(pi,len); // verkürzen und fertig
+	var cl_LF pires = sqrt(cl_I_to_LF(J3,actuallen)) / fsum;
+	return shorten(pires,len); // verkürzen und fertig
 }
 // Bit complexity (N := len): O(log(N)^2*M(N)).
 
@@ -256,22 +258,24 @@ const cl_LF compute_pi_ramanujan_163_fast (uintC len)
 //     it using binary splitting, is an O(log N * M(N)) algorithm, and
 //     outperforms all of the others.
 
-const cl_LF cl_pi (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);
 	if (len == oldlen)
 		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:
+	// > 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:
+	// 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