1 ///////////////////////////////////////////////////////////////
2 // Copyright 2013 John Maddock. Distributed under the Boost
3 // Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_
5
6 #ifndef BOOST_MULTIPRECISION_CPP_BIN_FLOAT_TRANSCENDENTAL_HPP
7 #define BOOST_MULTIPRECISION_CPP_BIN_FLOAT_TRANSCENDENTAL_HPP
8
9 namespace boost{ namespace multiprecision{ namespace backends{
10
11 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
eval_exp_taylor(cpp_bin_float<Digits,DigitBase,Allocator,Exponent,MinE,MaxE> & res,const cpp_bin_float<Digits,DigitBase,Allocator,Exponent,MinE,MaxE> & arg)12 void eval_exp_taylor(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
13 {
14 static const int bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count;
15 //
16 // Taylor series for small argument, note returns exp(x) - 1:
17 //
18 res = limb_type(0);
19 cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> num(arg), denom, t;
20 denom = limb_type(1);
21 eval_add(res, num);
22
23 for(unsigned k = 2; ; ++k)
24 {
25 eval_multiply(denom, k);
26 eval_multiply(num, arg);
27 eval_divide(t, num, denom);
28 eval_add(res, t);
29 if(eval_is_zero(t) || (res.exponent() - bits > t.exponent()))
30 break;
31 }
32 }
33
34 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
eval_exp(cpp_bin_float<Digits,DigitBase,Allocator,Exponent,MinE,MaxE> & res,const cpp_bin_float<Digits,DigitBase,Allocator,Exponent,MinE,MaxE> & arg)35 void eval_exp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
36 {
37 //
38 // This is based on MPFR's method, let:
39 //
40 // n = floor(x / ln(2))
41 //
42 // Then:
43 //
44 // r = x - n ln(2) : 0 <= r < ln(2)
45 //
46 // We can reduce r further by dividing by 2^k, with k ~ sqrt(n),
47 // so if:
48 //
49 // e0 = exp(r / 2^k) - 1
50 //
51 // With e0 evaluated by taylor series for small arguments, then:
52 //
53 // exp(x) = 2^n (1 + e0)^2^k
54 //
55 // Note that to preserve precision we actually square (1 + e0) k times, calculating
56 // the result less one each time, i.e.
57 //
58 // (1 + e0)^2 - 1 = e0^2 + 2e0
59 //
60 // Then add the final 1 at the end, given that e0 is small, this effectively wipes
61 // out the error in the last step.
62 //
63 using default_ops::eval_multiply;
64 using default_ops::eval_subtract;
65 using default_ops::eval_add;
66 using default_ops::eval_convert_to;
67 using default_ops::eval_increment;
68
69 int type = eval_fpclassify(arg);
70 bool isneg = eval_get_sign(arg) < 0;
71 if(type == (int)FP_NAN)
72 {
73 res = arg;
74 errno = EDOM;
75 return;
76 }
77 else if(type == (int)FP_INFINITE)
78 {
79 res = arg;
80 if(isneg)
81 res = limb_type(0u);
82 else
83 res = arg;
84 return;
85 }
86 else if(type == (int)FP_ZERO)
87 {
88 res = limb_type(1);
89 return;
90 }
91 cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> t, n;
92 if(isneg)
93 {
94 t = arg;
95 t.negate();
96 eval_exp(res, t);
97 t.swap(res);
98 res = limb_type(1);
99 eval_divide(res, t);
100 return;
101 }
102
103 eval_divide(n, arg, default_ops::get_constant_ln2<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >());
104 eval_floor(n, n);
105 eval_multiply(t, n, default_ops::get_constant_ln2<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >());
106 eval_subtract(t, arg);
107 t.negate();
108 if (t.compare(default_ops::get_constant_ln2<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >()) > 0)
109 {
110 // There are some rare cases where the multiply rounds down leaving a remainder > ln2
111 // See https://github.com/boostorg/multiprecision/issues/120
112 eval_increment(n);
113 t = limb_type(0);
114 }
115 if(eval_get_sign(t) < 0)
116 {
117 // There are some very rare cases where arg/ln2 is an integer, and the subsequent multiply
118 // rounds up, in that situation t ends up negative at this point which breaks our invariants below:
119 t = limb_type(0);
120 }
121
122 Exponent k, nn;
123 eval_convert_to(&nn, n);
124
125 if (nn == (std::numeric_limits<Exponent>::max)())
126 {
127 // The result will necessarily oveflow:
128 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
129 return;
130 }
131
132 BOOST_ASSERT(t.compare(default_ops::get_constant_ln2<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >()) < 0);
133
134 k = nn ? Exponent(1) << (msb(nn) / 2) : 0;
135 k = (std::min)(k, (Exponent)(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count / 4));
136 eval_ldexp(t, t, -k);
137
138 eval_exp_taylor(res, t);
139 //
140 // Square 1 + res k times:
141 //
142 for(Exponent s = 0; s < k; ++s)
143 {
144 t.swap(res);
145 eval_multiply(res, t, t);
146 eval_ldexp(t, t, 1);
147 eval_add(res, t);
148 }
149 eval_add(res, limb_type(1));
150 eval_ldexp(res, res, nn);
151 }
152
153 }}} // namespaces
154
155 #endif
156
157