src/float/transcendental/cl_LF_atan_recip.cc

   1 // cl_atan_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 #include "cl_alloca.h"
  17
  18 #undef floor
  19 #include <cmath>
  20 #define floor cln_floor
  21
  22 namespace cln {
  23
  24 // Method:
  25 // See examples/atan_recip.cc for a comparison of the algorithms.
  26 // Here we take algorithm 2d. It's the fastest throughout the range.
  27
  28 const cl_LF cl_atan_recip (cl_I m, uintC len)
  29 {
  30         var uintC actuallen = len + 1;
  31         var cl_I m2 = m*m+1;
  32         var uintL N = (uintL)(0.69314718*intDsize*actuallen/::log(double_approx(m2))) + 1;
  33         CL_ALLOCA_STACK;
  34         var cl_I* pv = (cl_I*) cl_alloca(N*sizeof(cl_I));
  35         var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
  36         var uintL n;
  37         new (&pv[0]) cl_I (m);
  38         new (&qv[0]) cl_I (m2);
  39         for (n = 1; n < N; n++) {
  40                 new (&pv[n]) cl_I (2*n);
  41                 new (&qv[n]) cl_I ((2*n+1)*m2);
  42         }
  43         var cl_pq_series series;
  44         series.pv = pv; series.qv = qv; series.qsv = NULL;
  45         var cl_LF result = eval_rational_series(N,series,actuallen);
  46         for (n = 0; n < N; n++) {
  47                 pv[n].~cl_I();
  48                 qv[n].~cl_I();
  49         }
  50         return shorten(result,len);
  51 }
  52 // Bit complexity (N = len): O(log(N)^2*M(N)).
  53
  54 }  // namespace cln