Replace unused macro with cl_unused.

[cln.git] / src / float / transcendental / cl_LF_exp1.cc

diff --git a/src/float/transcendental/cl_LF_exp1.cc b/src/float/transcendental/cl_LF_exp1.cc
index c9272fb434f444edf60538c2201f4e7a12553c72..c1f1acee5a98de5203bf0114a56e5f1fedbc079c 100644 (file)
--- a/src/float/transcendental/cl_LF_exp1.cc
+++ b/src/float/transcendental/cl_LF_exp1.cc
@@ -1,31 +1,32 @@
-// cl_exp1().
+// 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 "cl_lfloat.h"
-#include "cl_LF_tran.h"
-#include "cl_LF.h"
-#include "cl_integer.h"
-#include "cl_alloca.h"
+#include "cln/lfloat.h"
+#include "float/transcendental/cl_LF_tran.h"
+#include "float/lfloat/cl_LF.h"
+#include "cln/integer.h"
 
 #undef floor
-#include <math.h>
+#include <cmath>
 #define floor cln_floor
 
+namespace cln {
+
 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)
@@ -36,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)).
 
@@ -73,22 +76,24 @@ const cl_LF compute_exp1 (uintC len)
 // 25000   111
 // 50000   254
 
-const cl_LF cl_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