src/float/transcendental/cl_LF_atanh_recip.cc

   1 // cl_atanh_recip().
   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/integer.h"
  13 #include "cln/lfloat.h"
  14 #include "cl_LF.h"
  15 #include "cl_LF_tran.h"
  16
  17 #undef floor
  18 #include <cmath>
  19 #define floor cln_floor
  20
  21 namespace cln {
  22
  23 // Method:
  24 // See examples/atanh_recip.cc for a comparison of the algorithms.
  25 // Here we take algorithm 1d. It's the fastest throughout the range.
  26
  27 const cl_LF cl_atanh_recip (cl_I m, uintC len)
  28 {
  29         var uintC actuallen = len + 1;
  30         var uintC N = (uintC)(0.69314718*intDsize/2*actuallen/::log(double_approx(m))) + 1;
  31         struct rational_series_stream : cl_qb_series_stream {
  32                 var uintC n;
  33                 var cl_I m;
  34                 var cl_I m2;
  35                 static cl_qb_series_term computenext (cl_qb_series_stream& thisss)
  36                 {
  37                         var rational_series_stream& thiss = (rational_series_stream&)thisss;
  38                         var uintC n = thiss.n;
  39                         var cl_qb_series_term result;
  40                         result.b = 2*n+1;
  41                         result.q = (n==0 ? thiss.m : thiss.m2);
  42                         thiss.n = n+1;
  43                         return result;
  44                 }
  45                 rational_series_stream(const cl_I& m_)
  46                         : cl_qb_series_stream (rational_series_stream::computenext),
  47                           n(0), m(m_), m2(square(m_)) {}
  48         } series(m);
  49         var cl_LF result = eval_rational_series<false>(N,series,actuallen);
  50         return shorten(result,len);
  51 }
  52 // Bit complexity (N = len): O(log(N)^2*M(N)).
  53
  54 }  // namespace cln