]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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> | |
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> | |
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 | ||
68 | int type = eval_fpclassify(arg); | |
69 | bool isneg = eval_get_sign(arg) < 0; | |
70 | if(type == (int)FP_NAN) | |
71 | { | |
72 | res = arg; | |
73 | return; | |
74 | } | |
75 | else if(type == (int)FP_INFINITE) | |
76 | { | |
77 | res = arg; | |
78 | if(isneg) | |
79 | res = limb_type(0u); | |
80 | else | |
81 | res = arg; | |
82 | return; | |
83 | } | |
84 | else if(type == (int)FP_ZERO) | |
85 | { | |
86 | res = limb_type(1); | |
87 | return; | |
88 | } | |
89 | cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> t, n; | |
90 | if(isneg) | |
91 | { | |
92 | t = arg; | |
93 | t.negate(); | |
94 | eval_exp(res, t); | |
95 | t.swap(res); | |
96 | res = limb_type(1); | |
97 | eval_divide(res, t); | |
98 | return; | |
99 | } | |
100 | ||
101 | eval_divide(n, arg, default_ops::get_constant_ln2<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >()); | |
102 | eval_floor(n, n); | |
103 | eval_multiply(t, n, default_ops::get_constant_ln2<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >()); | |
104 | eval_subtract(t, arg); | |
105 | t.negate(); | |
106 | if(eval_get_sign(t) < 0) | |
107 | { | |
108 | // There are some very rare cases where arg/ln2 is an integer, and the subsequent multiply | |
109 | // rounds up, in that situation t ends up negative at this point which breaks our invariants below: | |
110 | t = limb_type(0); | |
111 | } | |
112 | BOOST_ASSERT(t.compare(default_ops::get_constant_ln2<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >()) < 0); | |
113 | ||
114 | Exponent k, nn; | |
115 | eval_convert_to(&nn, n); | |
116 | k = nn ? Exponent(1) << (msb(nn) / 2) : 0; | |
117 | eval_ldexp(t, t, -k); | |
118 | ||
119 | eval_exp_taylor(res, t); | |
120 | // | |
121 | // Square 1 + res k times: | |
122 | // | |
123 | for(int s = 0; s < k; ++s) | |
124 | { | |
125 | t.swap(res); | |
126 | eval_multiply(res, t, t); | |
127 | eval_ldexp(t, t, 1); | |
128 | eval_add(res, t); | |
129 | } | |
130 | eval_add(res, limb_type(1)); | |
131 | eval_ldexp(res, res, nn); | |
132 | } | |
133 | ||
134 | }}} // namespaces | |
135 | ||
136 | #endif | |
137 |