]> www.ginac.de Git - cln.git/blob - src/float/transcendental/cl_LF_atanh_recip.cc
Update to recently found large Mersenne prime.
[cln.git] / src / float / transcendental / cl_LF_atanh_recip.cc
1 // cl_atanh_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/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