X-Git-Url: https://ginac.de/CLN/cln.git//cln.git?a=blobdiff_plain;f=src%2Ffloat%2Ftranscendental%2Fcl_LF_exp1.cc;h=c1f1acee5a98de5203bf0114a56e5f1fedbc079c;hb=3af2cde18b3aabed4c808b0113daa81c2263b0bd;hp=713e2dde2c8f86e0f27b1605b13bbbaa6ac8563c;hpb=850abfde7f0d985ba01526c346bcd0d733562943;p=cln.git

diff --git a/src/float/transcendental/cl_LF_exp1.cc b/src/float/transcendental/cl_LF_exp1.cc
index 713e2dd..c1f1ace 100644
--- a/src/float/transcendental/cl_LF_exp1.cc
+++ b/src/float/transcendental/cl_LF_exp1.cc
@@ -1,19 +1,18 @@
 // exp1().
 
 // 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"
 
 #undef floor
 #include <cmath>
@@ -26,8 +25,8 @@ const cl_LF compute_exp1 (uintC len)
 	// Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
 	// with appropriate N, and
 	//   a(n) = 1, b(n) = 1, p(n) = 1, q(n) = n for n>0.
-	var uintC actuallen = len+1; // 1 Schutz-Digit
-	// How many terms to we need for M bits of precision? N terms suffice,
+	var uintC actuallen = len+1; // 1 guard digit
+	// How many terms do we need for M bits of precision? N terms suffice,
 	// provided that
 	//   1/N! < 2^-M
 	// <==   N*(log(N)-1) > M*log(2)
@@ -38,25 +37,27 @@ const cl_LF compute_exp1 (uintC len)
 	// Third approximation:
 	//   N2 = ceiling(M*log(2)/(log(N1)-1)), slightly too large.
 	//   N = N2+2, two more terms for safety.
-	var uintL N0 = intDsize*actuallen;
-	var uintL N1 = (uintL)(0.693147*intDsize*actuallen/(::log((double)N0)-1.0));
-	var uintL N2 = (uintL)(0.693148*intDsize*actuallen/(::log((double)N1)-1.0))+1;
-	var uintL N = N2+2;
-	CL_ALLOCA_STACK;
-	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, qv[n]) (n==0 ? 1 : n);
-	}
-	var cl_q_series series;
-	series.qv = qv;
-	series.qsv = (len >= 1000 && len <= 10000 ? qsv : 0); // tiny speedup
-	var cl_LF fsum = eval_rational_series(N,series,actuallen);
-	for (n = 0; n < N; n++) {
-		qv[n].~cl_I();
-	}
-	return shorten(fsum,len); // verkürzen und fertig
+	var uintC N0 = intDsize*actuallen;
+	var uintC N1 = (uintC)(0.693147*intDsize*actuallen/(::log((double)N0)-1.0));
+	var uintC N2 = (uintC)(0.693148*intDsize*actuallen/(::log((double)N1)-1.0))+1;
+	var uintC N = N2+2;
+	struct rational_series_stream : cl_q_series_stream {
+		var uintC n;
+		static cl_q_series_term computenext (cl_q_series_stream& thisss)
+		{
+			var rational_series_stream& thiss = (rational_series_stream&)thisss;
+			var uintC n = thiss.n;
+			var cl_q_series_term result;
+			result.q = (n==0 ? 1 : n);
+			thiss.n = n+1;
+			return result;
+		}
+		rational_series_stream()
+			: cl_q_series_stream (rational_series_stream::computenext),
+			  n(0) {}
+	} series;
+	var cl_LF fsum = eval_rational_series<false>(N,series,actuallen);
+	return shorten(fsum,len); // verkürzen und fertig
 }
 // Bit complexity (N = len): O(log(N)*M(N)).
 
@@ -77,22 +78,22 @@ const cl_LF compute_exp1 (uintC len)
 
 const cl_LF exp1 (uintC len)
 {
-	var uintC oldlen = TheLfloat(cl_LF_exp1)->len; // vorhandene Länge
+	var uintC oldlen = TheLfloat(cl_LF_exp1())->len; // vorhandene Länge
 	if (len < oldlen)
-		return shorten(cl_LF_exp1,len);
+		return shorten(cl_LF_exp1(),len);
 	if (len == oldlen)
-		return cl_LF_exp1;
+		return cl_LF_exp1();
 
-	// TheLfloat(cl_LF_exp1)->len um mindestens einen konstanten Faktor
-	// > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird:
+	// TheLfloat(cl_LF_exp1())->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_exp1 = compute_exp1(newlen); // (exp 1)
-	return (len < newlen ? shorten(cl_LF_exp1,len) : cl_LF_exp1);
+	// gewünschte > vorhandene Länge -> muß nachberechnen:
+	cl_LF_exp1() = compute_exp1(newlen); // (exp 1)
+	return (len < newlen ? shorten(cl_LF_exp1(),len) : cl_LF_exp1());
 }
 
 }  // namespace cln