src/float/transcendental/cl_LF_atan_recip.cc

   1 // cl_atan_recip().
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "float/transcendental/cl_F_tran.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cln/integer.h"
  13 #include "cln/lfloat.h"
  14 #include "float/lfloat/cl_LF.h"
  15 #include "float/transcendental/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/atan_recip.cc for a comparison of the algorithms.
  25 // Here we take algorithm 2d. It's the fastest throughout the range.
  26
  27 const cl_LF cl_atan_recip (cl_I m, uintC len)
  28 {
  29         var uintC actuallen = len + 1;
  30         var cl_I m2 = m*m+1;
  31         var uintC N = (uintC)(0.69314718*intDsize*actuallen/::log(double_approx(m2))) + 1;
  32         struct rational_series_stream : cl_pq_series_stream {
  33                 var uintC n;
  34                 var cl_I m;
  35                 var cl_I m2;
  36                 static cl_pq_series_term computenext (cl_pq_series_stream& thisss)
  37                 {
  38                         var rational_series_stream& thiss = (rational_series_stream&)thisss;
  39                         var uintC n = thiss.n;
  40                         var cl_pq_series_term result;
  41                         if (n==0) {
  42                                 result.p = thiss.m;
  43                                 result.q = thiss.m2;
  44                         } else {
  45                                 result.p = 2*n;
  46                                 result.q = (2*n+1)*thiss.m2;
  47                         }
  48                         thiss.n = n+1;
  49                         return result;
  50                 }
  51                 rational_series_stream(const cl_I& m_, const cl_I& m2_)
  52                         : cl_pq_series_stream (rational_series_stream::computenext),
  53                           n(0), m(m_), m2(m2_) {}
  54         } series(m,m2);
  55         var cl_LF result = eval_rational_series<false>(N,series,actuallen);
  56         return shorten(result,len);
  57 }
  58 // Bit complexity (N = len): O(log(N)^2*M(N)).
  59
  60 }  // namespace cln