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