More memory efficient Euler-Mascheroni constant:

* src/float/transcendental/cl_LF_tran.h (cl_pqd_series_stream): New.
* src/float/transcendental/cl_LF_ratsumseries_stream_pqd.cc: New file.
* src/float/transcendental/cl_LF_ratsumseries_stream_pqd_aux.cc: New file.
* src/float/transcendental/cl_LF_eulerconst.cc: Compute series coefficients
on demand, using a series stream object.

ChangeLog

@@ -1,3 +1,13 @@
+2007-09-07  Richard B. Kreckel  <kreckel@ginac.de>
+
+	More memory efficient Euler-Mascheroni constant:
+	* src/float/transcendental/cl_LF_tran.h (cl_pqd_series_stream): New.
+	* src/float/transcendental/cl_LF_ratsumseries_stream_pqd.cc: New file.
+	* src/float/transcendental/cl_LF_ratsumseries_stream_pqd_aux.cc: New 
+	file.
+	* src/float/transcendental/cl_LF_eulerconst.cc: Compute series
+	coefficients on demand, using a series stream object.
+
 2007-08-02  Richard B. Kreckel  <kreckel@ginac.de>
 
 	* src/base/digitseq/cl_DS_div.cc (cl_recip_suitable): uintC arguments.

src/float/transcendental/cl_LF_eulerconst.cc

@@ -447,16 +447,26 @@ const cl_LF compute_eulerconst_besselintegral4 (uintC len)
 	var uintC sx = (uintC)(0.25*0.693148*intDsize*actuallen)+1;
 	var uintC N = (uintC)(3.591121477*sx);
 	var cl_I x = square((cl_I)sx);
-	CL_ALLOCA_STACK;
-	var cl_pqd_series_term* args = (cl_pqd_series_term*) cl_alloca(N*sizeof(cl_pqd_series_term));
-	var uintC n;
-	for (n = 0; n < N; n++) {
-		init1(cl_I, args[n].p) (x);
-		init1(cl_I, args[n].q) (square((cl_I)(n+1)));
-		init1(cl_I, args[n].d) (n+1);
+	struct rational_series_stream : cl_pqd_series_stream {
+		uintC n;
+		cl_I x;
+		static cl_pqd_series_term computenext (cl_pqd_series_stream& thisss)
+		{
+			var rational_series_stream& thiss = (rational_series_stream&)thisss;
+			var uintC n = thiss.n;
+			var cl_pqd_series_term result;
+			result.p = thiss.x;
+			result.q = square((cl_I)(n+1));
+			result.d = n+1;
+			thiss.n = n+1;
+			return result;
 		}
+		rational_series_stream (const cl_I& _x)
+			: cl_pqd_series_stream (rational_series_stream::computenext),
+			  n (0), x (_x) {}
+	} series(x);
 	var cl_pqd_series_result sums;
-	eval_pqd_series_aux(N,args,sums);
+	eval_pqd_series_aux(N,series,sums);
 	// Instead of computing  fsum = 1 + T/Q  and  gsum = V/(D*Q)
 	// and then dividing them, to compute  gsum/fsum, we save two
 	// divisions by computing  V/(D*(Q+T)).
@@ -464,11 +474,6 @@ const cl_LF compute_eulerconst_besselintegral4 (uintC len)
 	  cl_I_to_LF(sums.V,actuallen)
 	  / The(cl_LF)(sums.D * cl_I_to_LF(sums.Q+sums.T,actuallen))
 	  - ln(cl_I_to_LF(sx,actuallen));
-	for (n = 0; n < N; n++) {
-		args[n].p.~cl_I();
-		args[n].q.~cl_I();
-		args[n].d.~cl_I();
-	}
 	return shorten(result,len); // verkürzen und fertig
 }
 // Bit complexity (N = len): O(log(N)^2*M(N)).

src/float/transcendental/cl_LF_ratsumseries_stream_pqd.cc

@@ -0,0 +1,31 @@
+// eval_pqd_series().
+
+// General includes.
+#include "cl_sysdep.h"
+
+// Specification.
+#include "cl_LF_tran.h"
+
+
+// Implementation.
+
+#include "cln/lfloat.h"
+#include "cln/integer.h"
+#include "cl_LF.h"
+
+namespace cln {
+
+const cl_LF eval_pqd_series (uintC N, cl_pqd_series_stream& args, uintC len)
+{
+	if (N==0)
+		return cl_I_to_LF(0,len);
+	var cl_pqd_series_result sums;
+	eval_pqd_series_aux(N,args,sums);
+	// Instead of computing  fsum = T/Q  and  gsum = V/(D*Q)
+	// and then dividing them, to compute  gsum/fsum, we save two
+	// divisions by computing  V/(D*T).
+	return
+	  cl_I_to_LF(sums.V,len) / The(cl_LF)(sums.D * cl_I_to_LF(sums.T,len));
+}
+
+}  // namespace cln

src/float/transcendental/cl_LF_ratsumseries_stream_pqd_aux.cc

@@ -0,0 +1,86 @@
+// eval_pqd_series_aux().
+
+// General includes.
+#include "cl_sysdep.h"
+
+// Specification.
+#include "cl_LF_tran.h"
+
+
+// Implementation.
+
+#include "cln/integer.h"
+#include "cln/exception.h"
+
+namespace cln {
+
+void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result& Z, bool rightmost)
+{
+	// N = N2-N1
+	switch (N) {
+	case 0:
+		throw runtime_exception(); break;
+	case 1: {
+		var cl_pqd_series_term v0 = args.next(); // [N1]
+		if (!rightmost) { Z.P = v0.p; }
+		Z.Q = v0.q;
+		Z.T = v0.p;
+		if (!rightmost) { Z.C = 1; }
+		Z.D = v0.d;
+		Z.V = v0.p;
+		break;
+		}
+	case 2: {
+		var cl_pqd_series_term v0 = args.next(); // [N1]
+		var cl_pqd_series_term v1 = args.next(); // [N1+1]
+		var cl_I p01 = v0.p * v1.p;
+		if (!rightmost) { Z.P = p01; }
+		Z.Q = v0.q * v1.q;
+		var cl_I p0q1 = v0.p * v1.q + p01;
+		Z.T = p0q1;
+		if (!rightmost) { Z.C = v1.d + v0.d; }
+		Z.D = v0.d * v1.d;
+		Z.V = v1.d * p0q1 + v0.d * p01;
+		break;
+		}
+	case 3: {
+		var cl_pqd_series_term v0 = args.next(); // [N1]
+		var cl_pqd_series_term v1 = args.next(); // [N1+1]
+		var cl_pqd_series_term v2 = args.next(); // [N1+2]
+		var cl_I p01 = v0.p * v1.p;
+		var cl_I p012 = p01 * v2.p;
+		if (!rightmost) { Z.P = p012; }
+		Z.Q = v0.q * v1.q * v2.q;
+		var cl_I p0q1 = v0.p * v1.q + p01;
+		Z.T = v2.q * p0q1 + p012;
+		var cl_I d01 = v0.d * v1.d;
+		if (!rightmost) { Z.C = (v1.d + v0.d) * v2.d + d01; }
+		Z.D = d01 * v2.d;
+		Z.V = v2.d * (v2.q * (v1.d * p0q1 + v0.d * p01) + (v1.d + v0.d) * p012) + d01 * p012;
+		break;
+		}
+	default: {
+		var uintC Nm = N/2; // midpoint
+		// Compute left part.
+		var cl_pqd_series_result L;
+		eval_pqd_series_aux(Nm,args,L,false);
+		// Compute right part.
+		var cl_pqd_series_result R;
+		eval_pqd_series_aux(N-Nm,args,R,rightmost);
+		// Put together partial results.
+		if (!rightmost) { Z.P = L.P * R.P; }
+		Z.Q = L.Q * R.Q;
+		// Z.S = L.S + L.P/L.Q*R.S;
+		var cl_I tmp = L.P * R.T;
+		Z.T = R.Q * L.T + tmp;
+		if (!rightmost) { Z.C = L.C * R.D + L.D * R.C; }
+		Z.D = L.D * R.D;
+		// Z.U = L.U + L.C/L.D * L.P/L.Q * R.S + L.P/L.Q * R.U;
+		// Z.V = R.D * R.Q * L.V + R.D * L.C * L.P * R.T + L.D * L.P * R.V;
+		Z.V = R.D * (R.Q * L.V + L.C * tmp) + L.D * L.P * R.V;
+		break;
+		}
+	}
+}
+
+}  // namespace cln

src/float/transcendental/cl_LF_tran.h

@@ -255,9 +255,17 @@ struct cl_pqd_series_result {
 	cl_I D;
 	cl_I V;
 };
+struct cl_pqd_series_stream {
+	cl_pqd_series_term (*nextfn)(cl_pqd_series_stream&);
+	cl_pqd_series_term next () { return nextfn(*this); }
+	// Constructor.
+	cl_pqd_series_stream( cl_pqd_series_term (*n)(cl_pqd_series_stream&)) : nextfn (n) {}
+};
 extern void eval_pqd_series_aux (uintC N, cl_pqd_series_term* args, cl_pqd_series_result& Z, bool rightmost = true);
+extern void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result& Z, bool rightmost = true);
 // Ditto, but returns U/S.
 extern const cl_LF eval_pqd_series (uintC N, cl_pqd_series_term* args, uintC len);
+extern const cl_LF eval_pqd_series (uintC N, cl_pqd_series_stream& args, uintC len);
 
 }  // namespace cln