src/float/transcendental/cl_LF_exp1.cc

   1 // exp1().
   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 #undef floor
  19 #include <cmath>
  20 #define floor cln_floor
  21
  22 namespace cln {
  23
  24 const cl_LF compute_exp1 (uintC len)
  25 {
  26         // Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
  27         // with appropriate N, and
  28         //   a(n) = 1, b(n) = 1, p(n) = 1, q(n) = n for n>0.
  29         var uintC actuallen = len+1; // 1 Schutz-Digit
  30         // How many terms to we need for M bits of precision? N terms suffice,
  31         // provided that
  32         //   1/N! < 2^-M
  33         // <==   N*(log(N)-1) > M*log(2)
  34         // First approximation:
  35         //   N0 = M will suffice, so put N<=N0.
  36         // Second approximation:
  37         //   N1 = floor(M*log(2)/(log(N0)-1)), slightly too small, so put N>=N1.
  38         // Third approximation:
  39         //   N2 = ceiling(M*log(2)/(log(N1)-1)), slightly too large.
  40         //   N = N2+2, two more terms for safety.
  41         var uintL N0 = intDsize*actuallen;
  42         var uintL N1 = (uintL)(0.693147*intDsize*actuallen/(::log((double)N0)-1.0));
  43         var uintL N2 = (uintL)(0.693148*intDsize*actuallen/(::log((double)N1)-1.0))+1;
  44         var uintL N = N2+2;
  45         CL_ALLOCA_STACK;
  46         var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
  47         var uintL* qsv = (uintL*) cl_alloca(N*sizeof(uintL));
  48         var uintL n;
  49         for (n = 0; n < N; n++) {
  50                 init1(cl_I, qv[n]) (n==0 ? 1 : n);
  51         }
  52         var cl_q_series series;
  53         series.qv = qv;
  54         series.qsv = (len >= 1000 && len <= 10000 ? qsv : 0); // tiny speedup
  55         var cl_LF fsum = eval_rational_series(N,series,actuallen);
  56         for (n = 0; n < N; n++) {
  57                 qv[n].~cl_I();
  58         }
  59         return shorten(fsum,len); // verkürzen und fertig
  60 }
  61 // Bit complexity (N = len): O(log(N)*M(N)).
  62
  63 // Timings of the above algorithm, on an i486 33 MHz, running Linux.
  64 //    N
  65 //    10     0.0040
  66 //    25     0.0096
  67 //    50     0.0218
  68 //   100     0.057
  69 //   250     0.24
  70 //   500     0.78
  71 //  1000     2.22
  72 //  2500     7.6
  73 //  5000    17.8
  74 // 10000    41.4
  75 // 25000   111
  76 // 50000   254
  77
  78 const cl_LF exp1 (uintC len)
  79 {
  80         var uintC oldlen = TheLfloat(cl_LF_exp1)->len; // vorhandene Länge
  81         if (len < oldlen)
  82                 return shorten(cl_LF_exp1,len);
  83         if (len == oldlen)
  84                 return cl_LF_exp1;
  85
  86         // TheLfloat(cl_LF_exp1)->len um mindestens einen konstanten Faktor
  87         // > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird:
  88         var uintC newlen = len;
  89         oldlen += floor(oldlen,2); // oldlen * 3/2
  90         if (newlen < oldlen)
  91                 newlen = oldlen;
  92
  93         // gewünschte > vorhandene Länge -> muß nachberechnen:
  94         cl_LF_exp1 = compute_exp1(newlen); // (exp 1)
  95         return (len < newlen ? shorten(cl_LF_exp1,len) : cl_LF_exp1);
  96 }
  97
  98 }  // namespace cln