src/float/transcendental/cl_LF_catalanconst.cc

   1 // cl_catalanconst().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_F_tran.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_lfloat.h"
  13 #include "cl_LF_tran.h"
  14 #include "cl_LF.h"
  15 #include "cl_integer.h"
  16 #include "cl_alloca.h"
  17
  18 const cl_LF compute_catalanconst_ramanujan (uintC len)
  19 {
  20         // [Jonathan M. Borwein, Peter B. Borwein: Pi and the AGM.
  21         //  Wiley 1987. Section 11.3, exercise 16 g, p. 386]
  22         // G = 3/8 * sum(n=0..infty, n!^2 / (2n+1)!*(2n+1))
  23         //     + pi/8 * log(2+sqrt(3)).
  24         // Every summand gives 0.6 new decimal digits in precision.
  25         // The sum is best evaluated using fixed-point arithmetic,
  26         // so that the precision is reduced for the later summands.
  27         var uintL actuallen = len + 2; // 2 Schutz-Digits
  28         var sintL scale = intDsize*actuallen;
  29         var cl_I sum = 0;
  30         var cl_I n = 0;
  31         var cl_I factor = ash(1,scale);
  32         while (!zerop(factor)) {
  33                 sum = sum + truncate1(factor,2*n+1);
  34                 n = n+1;
  35                 factor = truncate1(factor*n,2*(2*n+1));
  36         }
  37         var cl_LF fsum = scale_float(cl_I_to_LF(sum,actuallen),-scale);
  38         var cl_LF g =
  39           scale_float(The(cl_LF)(3*fsum)
  40                       + The(cl_LF)(cl_pi(actuallen))
  41                         * ln(cl_I_to_LF(2,actuallen)+sqrt(cl_I_to_LF(3,actuallen))),
  42                       -3);
  43         return shorten(g,len); // verkürzen und fertig
  44 }
  45 // Bit complexity (N := len): O(N^2).
  46
  47 const cl_LF compute_catalanconst_ramanujan_fast (uintC len)
  48 {
  49         // Same formula as above, using a binary splitting evaluation.
  50         // See [Borwein, Borwein, section 10.2.3].
  51         var uintL actuallen = len + 2; // 2 Schutz-Digits
  52         // Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
  53         // with appropriate N, and
  54         //   a(n) = 1, b(n) = 2*n+1,
  55         //   p(n) = n for n>0, q(n) = 2*(2*n+1) for n>0.
  56         var uintL N = (intDsize/2)*actuallen;
  57         // 4^-N <= 2^(-intDsize*actuallen).
  58         CL_ALLOCA_STACK;
  59         var cl_I* bv = (cl_I*) cl_alloca(N*sizeof(cl_I));
  60         var cl_I* pv = (cl_I*) cl_alloca(N*sizeof(cl_I));
  61         var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
  62         var uintL n;
  63         init1(cl_I, bv[0]) (1);
  64         init1(cl_I, pv[0]) (1);
  65         init1(cl_I, qv[0]) (1);
  66         for (n = 1; n < N; n++) {
  67                 init1(cl_I, bv[n]) (2*n+1);
  68                 init1(cl_I, pv[n]) (n);
  69                 init1(cl_I, qv[n]) (2*(2*n+1));
  70         }
  71         var cl_pqb_series series;
  72         series.bv = bv;
  73         series.pv = pv; series.qv = qv; series.qsv = NULL;
  74         var cl_LF fsum = eval_rational_series(N,series,actuallen);
  75         for (n = 0; n < N; n++) {
  76                 bv[n].~cl_I();
  77                 pv[n].~cl_I();
  78                 qv[n].~cl_I();
  79         }
  80         var cl_LF g =
  81           scale_float(The(cl_LF)(3*fsum)
  82                       + The(cl_LF)(cl_pi(actuallen))
  83                         * ln(cl_I_to_LF(2,actuallen)+sqrt(cl_I_to_LF(3,actuallen))),
  84                       -3);
  85         return shorten(g,len); // verkürzen und fertig
  86 }
  87 // Bit complexity (N := len): O(log(N)^2*M(N)).
  88
  89 const cl_LF compute_catalanconst_expintegral1 (uintC len)
  90 {
  91         // We compute f(x) classically and g(x) using the partial sums of f(x).
  92         var uintC actuallen = len+2; // 2 Schutz-Digits
  93         var uintL x = (uintL)(0.693148*intDsize*actuallen)+1;
  94         var uintL N = (uintL)(2.718281828*x);
  95         var cl_LF fterm = cl_I_to_LF(1,actuallen);
  96         var cl_LF fsum = fterm;
  97         var cl_LF gterm = fterm;
  98         var cl_LF gsum = gterm;
  99         var uintL n;
 100         // After n loops
 101         //   fterm = x^n/n!, fsum = 1 + x/1! + ... + x^n/n!,
 102         //   gterm = S_n*x^n/n!, gsum = S_0*x^0/0! + ... + S_n*x^n/n!.
 103         for (n = 1; n < N; n++) {
 104                 fterm = The(cl_LF)(fterm*x)/n;
 105                 fsum = fsum + fterm;
 106                 gterm = The(cl_LF)(gterm*x)/n;
 107                 if (evenp(n))
 108                         gterm = gterm + fterm/square((cl_I)(2*n+1));
 109                 else
 110                         gterm = gterm - fterm/square((cl_I)(2*n+1));
 111                 gsum = gsum + gterm;
 112         }
 113         var cl_LF result = gsum/fsum;
 114         return shorten(result,len); // verkürzen und fertig
 115 }
 116 // Bit complexity (N = len): O(N^2).
 117
 118 // Same algorithm as expintegral1, but using binary splitting to evaluate
 119 // the sums.
 120 const cl_LF compute_catalanconst_expintegral2 (uintC len)
 121 {
 122         var uintC actuallen = len+2; // 2 Schutz-Digits
 123         var uintL x = (uintL)(0.693148*intDsize*actuallen)+1;
 124         var uintL N = (uintL)(2.718281828*x);
 125         CL_ALLOCA_STACK;
 126         var cl_pqd_series_term* args = (cl_pqd_series_term*) cl_alloca(N*sizeof(cl_pqd_series_term));
 127         var uintL n;
 128         for (n = 0; n < N; n++) {
 129                 if (n==0) {
 130                         init1(cl_I, args[n].p) (1);
 131                         init1(cl_I, args[n].q) (1);
 132                 } else {
 133                         init1(cl_I, args[n].p) (x);
 134                         init1(cl_I, args[n].q) (n);
 135                 }
 136                 init1(cl_I, args[n].d) (evenp(n)
 137                                         ? square((cl_I)(2*n+1))
 138                                         : -square((cl_I)(2*n+1)));
 139         }
 140         var cl_LF result = eval_pqd_series(N,args,actuallen);
 141         for (n = 0; n < N; n++) {
 142                 args[n].p.~cl_I();
 143                 args[n].q.~cl_I();
 144                 args[n].d.~cl_I();
 145         }
 146         return shorten(result,len); // verkürzen und fertig
 147 }
 148 // Bit complexity (N = len): O(log(N)^2*M(N)).
 149
 150 // Using Cohen-Villegas-Zagier acceleration, but without binary splitting.
 151 const cl_LF compute_catalanconst_cvz1 (uintC len)
 152 {
 153         var uintC actuallen = len+2; // 2 Schutz-Digits
 154         var uintL N = (uintL)(0.39321985*intDsize*actuallen)+1;
 155 #if 0
 156         var cl_LF fterm = cl_I_to_LF(2*(cl_I)N*(cl_I)N,actuallen);
 157         var cl_LF fsum = fterm;
 158         var cl_LF gterm = fterm;
 159         var cl_LF gsum = gterm;
 160         var uintL n;
 161         // After n loops
 162         //   fterm = (N+n)!N/(2n+2)!(N-n-1)!*2^(2n+2), fsum = ... + fterm,
 163         //   gterm = S_n*fterm, gsum = ... + gterm.
 164         for (n = 1; n < N; n++) {
 165                 fterm = The(cl_LF)(fterm*(2*(cl_I)(N-n)*(cl_I)(N+n)))/((cl_I)(2*n+1)*(cl_I)(n+1));
 166                 fsum = fsum + fterm;
 167                 gterm = The(cl_LF)(gterm*(2*(cl_I)(N-n)*(cl_I)(N+n)))/((cl_I)(2*n+1)*(cl_I)(n+1));
 168                 if (evenp(n))
 169                         gterm = gterm + fterm/square((cl_I)(2*n+1));
 170                 else
 171                         gterm = gterm - fterm/square((cl_I)(2*n+1));
 172                 gsum = gsum + gterm;
 173         }
 174         var cl_LF result = gsum/(cl_I_to_LF(1,actuallen)+fsum);
 175 #else
 176         // Take advantage of the fact that fterm and fsum are integers.
 177         var cl_I fterm = 2*(cl_I)N*(cl_I)N;
 178         var cl_I fsum = fterm;
 179         var cl_LF gterm = cl_I_to_LF(fterm,actuallen);
 180         var cl_LF gsum = gterm;
 181         var uintL n;
 182         // After n loops
 183         //   fterm = (N+n)!N/(2n+2)!(N-n-1)!*2^(2n+2), fsum = ... + fterm,
 184         //   gterm = S_n*fterm, gsum = ... + gterm.
 185         for (n = 1; n < N; n++) {
 186                 fterm = exquopos(fterm*(2*(cl_I)(N-n)*(cl_I)(N+n)),(cl_I)(2*n+1)*(cl_I)(n+1));
 187                 fsum = fsum + fterm;
 188                 gterm = The(cl_LF)(gterm*(2*(cl_I)(N-n)*(cl_I)(N+n)))/((cl_I)(2*n+1)*(cl_I)(n+1));
 189                 if (evenp(n))
 190                         gterm = gterm + cl_I_to_LF(fterm,actuallen)/square((cl_I)(2*n+1));
 191                 else
 192                         gterm = gterm - cl_I_to_LF(fterm,actuallen)/square((cl_I)(2*n+1));
 193                 gsum = gsum + gterm;
 194         }
 195         var cl_LF result = gsum/cl_I_to_LF(1+fsum,actuallen);
 196 #endif
 197         return shorten(result,len); // verkürzen und fertig
 198 }
 199 // Bit complexity (N = len): O(N^2).
 200
 201 // Using Cohen-Villegas-Zagier acceleration, with binary splitting.
 202 const cl_LF compute_catalanconst_cvz2 (uintC len)
 203 {
 204         var uintC actuallen = len+2; // 2 Schutz-Digits
 205         var uintL N = (uintL)(0.39321985*intDsize*actuallen)+1;
 206         CL_ALLOCA_STACK;
 207         var cl_pqd_series_term* args = (cl_pqd_series_term*) cl_alloca(N*sizeof(cl_pqd_series_term));
 208         var uintL n;
 209         for (n = 0; n < N; n++) {
 210                 init1(cl_I, args[n].p) (2*(cl_I)(N-n)*(cl_I)(N+n));
 211                 init1(cl_I, args[n].q) ((cl_I)(2*n+1)*(cl_I)(n+1));
 212                 init1(cl_I, args[n].d) (evenp(n)
 213                                         ? square((cl_I)(2*n+1))
 214                                         : -square((cl_I)(2*n+1)));
 215         }
 216         var cl_pqd_series_result sums;
 217         eval_pqd_series_aux(N,args,sums);
 218         // Here we need U/(1+S) = V/D(Q+T).
 219         var cl_LF result =
 220           cl_I_to_LF(sums.V,actuallen) / The(cl_LF)(sums.D * cl_I_to_LF(sums.Q+sums.T,actuallen));
 221         for (n = 0; n < N; n++) {
 222                 args[n].p.~cl_I();
 223                 args[n].q.~cl_I();
 224                 args[n].d.~cl_I();
 225         }
 226         return shorten(result,len); // verkürzen und fertig
 227 }
 228 // Bit complexity (N = len): O(log(N)^2*M(N)).
 229
 230 // Timings of the above algorithms, on an i486 33 MHz, running Linux.
 231 //    N      ram    ramfast  exp1    exp2    cvz1    cvz2
 232 //    10     0.055   0.068   0.32    0.91    0.076   0.11
 233 //    25     0.17    0.26    0.95    3.78    0.23    0.43
 234 //    50     0.43    0.73    2.81   11.5     0.63    1.36
 235 //   100     1.32    2.24    8.82   34.1     1.90    4.48
 236 //   250     6.60   10.4    48.7   127.5    10.3    20.8
 237 //   500    24.0    30.9   186     329      38.4    58.6
 238 //  1000    83.0    89.0   944     860     149     163
 239 //  2500   446     352    6964    3096    1032     545
 240 //  5000  1547     899
 241 // asymp.    N^2     FAST    N^2     FAST    N^2     FAST
 242 // (FAST means O(log(N)^2*M(N)))
 243 //
 244 // The "exp1" and "exp2" algorithms are always about 10 to 15 times slower
 245 // than the "ram" and "ramfast" algorithms.
 246 // The break-even point between "ram" and "ramfast" is at about N = 1410.
 247
 248 const cl_LF compute_catalanconst (uintC len)
 249 {
 250         if (len >= 1410)
 251                 return compute_catalanconst_ramanujan_fast(len);
 252         else
 253                 return compute_catalanconst_ramanujan(len);
 254 }
 255 // Bit complexity (N := len): O(log(N)^2*M(N)).
 256
 257 const cl_LF cl_catalanconst (uintC len)
 258 {
 259         var uintC oldlen = TheLfloat(cl_LF_catalanconst)->len; // vorhandene Länge
 260         if (len < oldlen)
 261                 return shorten(cl_LF_catalanconst,len);
 262         if (len == oldlen)
 263                 return cl_LF_catalanconst;
 264
 265         // TheLfloat(cl_LF_catalanconst)->len um mindestens einen konstanten Faktor
 266         // > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird:
 267         var uintC newlen = len;
 268         oldlen += floor(oldlen,2); // oldlen * 3/2
 269         if (newlen < oldlen)
 270                 newlen = oldlen;
 271
 272         // gewünschte > vorhandene Länge -> muß nachberechnen:
 273         cl_LF_catalanconst = compute_catalanconst(newlen);
 274         return (len < newlen ? shorten(cl_LF_catalanconst,len) : cl_LF_catalanconst);
 275 }