Replace unused macro with cl_unused.

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

diff --git a/src/float/transcendental/cl_LF_zeta3.cc b/src/float/transcendental/cl_LF_zeta3.cc
index 0455da5e36b28df06dd26419f01c85c849d59edb..4e56cfe8d04e1e0ed3930972af3a49e8a40bf999 100644 (file)
--- a/src/float/transcendental/cl_LF_zeta3.cc
+++ b/src/float/transcendental/cl_LF_zeta3.cc
@@ -1,24 +1,47 @@
-// cl_zeta3().
+// zeta3().
 
 // 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"
+#include "base/cl_alloca.h"
 
-const cl_LF cl_zeta3 (uintC len)
+namespace cln {
+
+const cl_LF zeta3 (uintC len)
 {
+       struct rational_series_stream : cl_pqa_series_stream {
+               uintC n;
+               static cl_pqa_series_term computenext (cl_pqa_series_stream& thisss)
+               {
+                       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;
+                       } else {
+                               result.p = -expt_pos(n,5);
+                       }
+                       result.q = expt_pos(2*n+1,5)<<5;
+                       result.a = 205*square((cl_I)n) + 250*(cl_I)n + 77;
+                       thiss.n = n+1;
+                       return result;
+               }
+               rational_series_stream ()
+                       : cl_pqa_series_stream (rational_series_stream::computenext),
+                         n (0) {}
+       } series;
        // Method:
-       //            /infinity                                  \ 
+       //            /infinity                                  \ 
        //            | -----       (n + 1)       2              |
        //        1   |  \      (-1)        (205 n  - 160 n + 32)|
        //        -   |   )     ---------------------------------|
@@ -27,40 +50,19 @@ const cl_LF cl_zeta3 (uintC len)
        //            \ n = 1                                    /
        //
        // The formula used to compute Zeta(3) has reference in the paper
-       // "Acceleration of Hypergeometric Series via the WZ method" by
-       // T. Amdeberhan and Doron Zeilberger, to appear in the Electronic
-       // Journal of Combinatorics [Wilf Festschrift Volume].
+       // "Hypergeometric Series Acceleration via the WZ method" by
+       // T. Amdeberhan and Doron Zeilberger,
+       // Electronic J. Combin. 4 (1997), R3.
        //
        // Computation of the sum:
        // Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
        // with appropriate N, and
        //   a(n) = 205*n^2+250*n+77, b(n) = 1,
        //   p(0) = 1, p(n) = -n^5 for n>0, q(n) = 32*(2n+1)^5.
-       var uintC actuallen = len+2; // 2 Schutz-Digits
-       var uintL N = ceiling(actuallen*intDsize,10);
+       var uintC actuallen = len+2; // 2 guard digits
+       var uintC N = ceiling(actuallen*intDsize,10);
        // 1024^-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 n;
-       for (n = 0; n < N; n++) {
-               init1(cl_I, av[n]) (205*square((cl_I)n) + 250*(cl_I)n + 77);
-               if (n==0)
-                       init1(cl_I, pv[n]) (1);
-               else
-                       init1(cl_I, pv[n]) (-expt_pos(n,5));
-               init1(cl_I, qv[n]) (expt_pos(2*n+1,5)<<5);
-       }
-       var cl_pqa_series series;
-       series.av = av;
-       series.pv = pv; series.qv = qv; series.qsv = NULL;
-       var cl_LF sum = 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 sum = eval_rational_series<false>(N,series,actuallen,actuallen);
        return scale_float(shorten(sum,len),-1);
 }
 // Bit complexity (N := len): O(log(N)^2*M(N)).
@@ -81,3 +83,5 @@ const cl_LF cl_zeta3 (uintC len)
 // 50000                          3723
 // asymp.    FAST     N^2    FAST    FAST
 // (FAST means O(log(N)^2*M(N)))
+
+}  // namespace cln