]> www.ginac.de Git - cln.git/blob - src/float/transcendental/cl_LF_coshsinh.cc
* include/cln/modules.h (CL_JUMP_TO): Fix AMD64 brokenness.
[cln.git] / src / float / transcendental / cl_LF_coshsinh.cc
1 // cl_coshsinh_ratseries().
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.h"
14 #include "cln/integer.h"
15
16 namespace cln {
17
18 inline const cl_LF_cosh_sinh_t operator* (const cl_LF_cosh_sinh_t& a, const cl_LF_cosh_sinh_t& b)
19 {
20         return cl_LF_cosh_sinh_t(a.cosh*b.cosh+a.sinh*b.sinh,a.sinh*b.cosh+a.cosh*b.sinh);
21 }
22
23 const cl_LF_cosh_sinh_t cl_coshsinh_ratseries (const cl_LF& x)
24 {
25         // Similar to expx_ratseries.
26         var uintC len = TheLfloat(x)->len;
27         var cl_idecoded_float x_ = integer_decode_float(x);
28         // x = (-1)^sign * 2^exponent * mantissa
29         var uintL lq = cl_I_to_UL(- x_.exponent);
30         var const cl_I& p = x_.mantissa;
31         // Compute sinh(p/2^lq) and cosh(p/2^lq) by splitting into pieces.
32         var cl_boolean first_factor = cl_true;
33         var cl_LF_cosh_sinh_t product;
34         var uintL b1;
35         var uintL b2;
36         for (b1 = 0, b2 = 1; b1 < lq; b1 = b2, b2 = 2*b2) {
37                 // Piece containing bits b1+1..b2 after "decimal point"
38                 // in the binary representation of (p/2^lq).
39                 var uintL lqk = (lq >= b2 ? b2 : lq);
40                 var cl_I pk = ldb(p,cl_byte(lqk-b1,lq-lqk));
41                 // Compute sinh(pk/2^lqk) and cosh(pk/2^lqk).
42                 if (!zerop(pk)) {
43                         if (minusp(x_.sign)) { pk = -pk; }
44                         var cl_LF_cosh_sinh_t factor = cl_coshsinh_aux(pk,lqk,len);
45                         if (first_factor) {
46                                 product = factor;
47                                 first_factor = cl_false;
48                         } else
49                                 product = product * factor;
50                 }
51         }
52         if (first_factor)
53                 return cl_LF_cosh_sinh_t(cl_I_to_LF(1,len),cl_I_to_LF(0,len));
54         else
55                 return product;
56 }
57 // Bit complexity (N = length(x)): O(log(N)^2*M(N)).
58
59 }  // namespace cln