src/float/transcendental/cl_LF_exp1.cc

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