src/float/transcendental/cl_LF_coshsinh_aux.cc

   1 // cl_coshsinh_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 "cl_alloca.h"
  17 #include "cln/abort.h"
  18
  19 #undef floor
  20 #include <cmath>
  21 #define floor cln_floor
  22
  23 namespace cln {
  24
  25 // Computing cosh(x) = sqrt(1+sinh(x)^2) instead of computing separately
  26 // by a power series evaluation brings 20% speedup, even more for small lengths.
  27 #define TRIVIAL_SPEEDUP
  28
  29 const cl_LF_cosh_sinh_t cl_coshsinh_aux (const cl_I& p, uintL lq, uintC len)
  30 {
  31  {      Mutable(cl_I,p);
  32         var uintL lp = integer_length(p); // now |p| < 2^lp.
  33         if (!(lp <= lq)) cl_abort();
  34         lp = lq - lp; // now |p/2^lq| < 2^-lp.
  35         // Minimize lq (saves computation time).
  36         {
  37                 var uintL lp2 = ord2(p);
  38                 if (lp2 > 0) {
  39                         p = p >> lp2;
  40                         lq = lq - lp2;
  41                 }
  42         }
  43         // cosh(p/2^lq):
  44         // Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
  45         // with appropriate N, and
  46         //   a(n) = 1, b(n) = 1,
  47         //   p(n) = p^2 for n>0,
  48         //   q(n) = (2*n-1)*(2*n)*(2^lq)^2 for n>0.
  49         // sinh(p/2^lq):
  50         // Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
  51         // with appropriate N, and
  52         //   a(n) = 1, b(n) = 1,
  53         //   p(0) = p, p(n) = p^2 for n>0,
  54         //   q(0) = 2^lq, q(n) = (2*n)*(2*n+1)*(2^lq)^2 for n>0.
  55         var uintC actuallen = len+1; // 1 Schutz-Digit
  56         // How many terms to we need for M bits of precision? N/2 terms suffice,
  57         // provided that
  58         //   1/(2^(N*lp)*N!) < 2^-M
  59         // <==   N*(log(N)-1)+N*lp*log(2) > M*log(2)
  60         // First approximation:
  61         //   N0 = M will suffice, so put N<=N0.
  62         // Second approximation:
  63         //   N1 = floor(M*log(2)/(log(N0)-1+lp*log(2))), slightly too small,
  64         //   so put N>=N1.
  65         // Third approximation:
  66         //   N2 = ceiling(M*log(2)/(log(N1)-1+lp*log(2))), slightly too large.
  67         //   N = N2+2, two more terms for safety.
  68         var uintL N0 = intDsize*actuallen;
  69         var uintL N1 = (uintL)(0.693147*intDsize*actuallen/(::log((double)N0)-1.0+0.693148*lp));
  70         var uintL N2 = (uintL)(0.693148*intDsize*actuallen/(::log((double)N1)-1.0+0.693147*lp))+1;
  71         var uintL N = N2+2;
  72         N = ceiling(N,2);
  73         CL_ALLOCA_STACK;
  74         var cl_I* pv = (cl_I*) cl_alloca(N*sizeof(cl_I));
  75         var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
  76         var uintL* qsv = (uintL*) cl_alloca(N*sizeof(uintL));
  77         var uintL n;
  78         var cl_I p2 = square(p);
  79         var cl_LF sinhsum;
  80         {
  81                 init1(cl_I, pv[0]) (p);
  82                 init1(cl_I, qv[0]) ((cl_I)1 << lq);
  83                 for (n = 1; n < N; n++) {
  84                         init1(cl_I, pv[n]) (p2);
  85                         init1(cl_I, qv[n]) (((cl_I)n*(cl_I)(2*n+1)) << (2*lq+1));
  86                 }
  87                 var cl_pq_series series;
  88                 series.pv = pv; series.qv = qv; series.qsv = qsv;
  89                 sinhsum = eval_rational_series(N,series,actuallen);
  90                 for (n = 0; n < N; n++) {
  91                         pv[n].~cl_I();
  92                         qv[n].~cl_I();
  93                 }
  94         }
  95         #if !defined(TRIVIAL_SPEEDUP)
  96         var cl_LF coshsum;
  97         {
  98                 init1(cl_I, pv[0]) (1);
  99                 init1(cl_I, qv[0]) (1);
 100                 for (n = 1; n < N; n++) {
 101                         init1(cl_I, pv[n]) (p2);
 102                         init1(cl_I, qv[n]) (((cl_I)n*(cl_I)(2*n-1)) << (2*lq+1));
 103                 }
 104                 var cl_pq_series series;
 105                 series.pv = pv; series.qv = qv; series.qsv = qsv;
 106                 coshsum = eval_rational_series(N,series,actuallen);
 107                 for (n = 0; n < N; n++) {
 108                         pv[n].~cl_I();
 109                         qv[n].~cl_I();
 110                 }
 111         }
 112         #else // TRIVIAL_SPEEDUP
 113         var cl_LF coshsum = sqrt(cl_I_to_LF(1,actuallen) + square(sinhsum));
 114         #endif
 115         return cl_LF_cosh_sinh_t(shorten(coshsum,len),shorten(sinhsum,len)); // verkürzen und fertig
 116 }}
 117 // Bit complexity (N = len, and if p has length O(log N) and ql = O(log N)):
 118 // O(log(N)*M(N)).
 119
 120 }  // namespace cln