src/float/transcendental/cl_LF_exp_aux.cc

   1 // cl_exp_aux().
   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 "cln/exception.h"
  17
  18 #undef floor
  19 #include <cmath>
  20 #define floor cln_floor
  21
  22 namespace cln {
  23
  24 const cl_LF cl_exp_aux (const cl_I& p, uintE lq, uintC len)
  25 {
  26  {      Mutable(cl_I,p);
  27         var uintE lp = integer_length(p); // now |p| < 2^lp.
  28         if (!(lp <= lq)) throw runtime_exception();
  29         lp = lq - lp; // now |p/2^lq| < 2^-lp.
  30         // Minimize lq (saves computation time).
  31         {
  32                 var uintC lp2 = ord2(p);
  33                 if (lp2 > 0) {
  34                         p = p >> lp2;
  35                         lq = lq - lp2;
  36                 }
  37         }
  38         // Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
  39         // with appropriate N, and
  40         //   a(n) = 1, b(n) = 1, p(n) = p for n>0, q(n) = n*2^lq for n>0.
  41         var uintC actuallen = len+1; // 1 guard digit
  42         // How many terms do we need for M bits of precision? N terms suffice,
  43         // provided that
  44         //   1/(2^(N*lp)*N!) < 2^-M
  45         // <==   N*(log(N)-1)+N*lp*log(2) > M*log(2)
  46         // First approximation:
  47         //   N0 = M will suffice, so put N<=N0.
  48         // Second approximation:
  49         //   N1 = floor(M*log(2)/(log(N0)-1+lp*log(2))), slightly too small,
  50         //   so put N>=N1.
  51         // Third approximation:
  52         //   N2 = ceiling(M*log(2)/(log(N1)-1+lp*log(2))), slightly too large.
  53         //   N = N2+2, two more terms for safety.
  54         var uintC N0 = intDsize*actuallen;
  55         var uintC N1 = (uintC)(0.693147*intDsize*actuallen/(::log((double)N0)-1.0+0.693148*lp));
  56         var uintC N2 = (uintC)(0.693148*intDsize*actuallen/(::log((double)N1)-1.0+0.693147*lp))+1;
  57         var uintC N = N2+2;
  58         struct rational_series_stream : cl_pq_series_stream {
  59                 var uintC n;
  60                 var cl_I p;
  61                 var uintE lq;
  62                 static cl_pq_series_term computenext (cl_pq_series_stream& thisss)
  63                 {
  64                         var rational_series_stream& thiss = (rational_series_stream&)thisss;
  65                         var uintC n = thiss.n;
  66                         var cl_pq_series_term result;
  67                         if (n==0) {
  68                                 result.p = 1;
  69                                 result.q = 1;
  70                         } else {
  71                                 result.p = thiss.p;
  72                                 result.q = (cl_I)n << thiss.lq;
  73                         }
  74                         thiss.n = n+1;
  75                         return result;
  76                 }
  77                 rational_series_stream(const cl_I& p_, uintE lq_)
  78                         : cl_pq_series_stream (rational_series_stream::computenext),
  79                           n (0), p(p_), lq(lq_) {}
  80         } series(p, lq);
  81         var cl_LF fsum = eval_rational_series<true>(N,series,actuallen);
  82         return shorten(fsum,len); // verkürzen und fertig
  83 }}
  84 // Bit complexity (N = len, and if p has length O(log N) and ql = O(log N)):
  85 // O(log(N)*M(N)).
  86
  87 }  // namespace cln