src/float/transcendental/cl_LF_catalanconst.cc

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