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

source code of boost/libs/multiprecision/include/boost/multiprecision/cpp_bin_float/transcendental.hpp