multiprecision/cpp_bin_float/transcendental.hpp

0001 ///////////////////////////////////////////////////////////////
0002 //  Copyright 2013 John Maddock. Distributed under the Boost
0003 //  Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
0005
0006 #ifndef BOOST_MP_CPP_BIN_FLOAT_TRANSCENDENTAL_HPP
0007 #define BOOST_MP_CPP_BIN_FLOAT_TRANSCENDENTAL_HPP
0008
0009 #include <boost/multiprecision/detail/assert.hpp>
0010
0011 namespace boost { namespace multiprecision { namespace backends {
0012
0013 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
0014 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)
0015 {
0016    constexpr std::ptrdiff_t bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count;
0017    //
0018    // Taylor series for small argument, note returns exp(x) - 1:
0019    //
0020    res = limb_type(0);
0021    cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> num(arg), denom, t;
0022    denom = limb_type(1);
0023    eval_add(res, num);
0024
0025    for (std::size_t k = 2;; ++k)
0026    {
0027       eval_multiply(denom, k);
0028       eval_multiply(num, arg);
0029       eval_divide(t, num, denom);
0030       eval_add(res, t);
0031       if (eval_is_zero(t) || (res.exponent() - bits > t.exponent()))
0032          break;
0033    }
0034 }
0035
0036 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
0037 void eval_exp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
0038 {
0039    //
0040    // This is based on MPFR's method, let:
0041    //
0042    // n = floor(x / ln(2))
0043    //
0044    // Then:
0045    //
0046    // r = x - n ln(2) : 0 <= r < ln(2)
0047    //
0048    // We can reduce r further by dividing by 2^k, with k ~ sqrt(n),
0049    // so if:
0050    //
0051    // e0 = exp(r / 2^k) - 1
0052    //
0053    // With e0 evaluated by taylor series for small arguments, then:
0054    //
0055    // exp(x) = 2^n (1 + e0)^2^k
0056    //
0057    // Note that to preserve precision we actually square (1 + e0) k times, calculating
0058    // the result less one each time, i.e.
0059    //
0060    // (1 + e0)^2 - 1 = e0^2 + 2e0
0061    //
0062    // Then add the final 1 at the end, given that e0 is small, this effectively wipes
0063    // out the error in the last step.
0064    //
0065    using default_ops::eval_add;
0066    using default_ops::eval_convert_to;
0067    using default_ops::eval_increment;
0068    using default_ops::eval_multiply;
0069    using default_ops::eval_subtract;
0070
0071    int  type  = eval_fpclassify(arg);
0072    bool isneg = eval_get_sign(arg) < 0;
0073    if (type == static_cast<int>(FP_NAN))
0074    {
0075       res   = arg;
0076       errno = EDOM;
0077       return;
0078    }
0079    else if (type == static_cast<int>(FP_INFINITE))
0080    {
0081       res = arg;
0082       if (isneg)
0083          res = limb_type(0u);
0084       else
0085          res = arg;
0086       return;
0087    }
0088    else if (type == static_cast<int>(FP_ZERO))
0089    {
0090       res = limb_type(1);
0091       return;
0092    }
0093    cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> t, n;
0094    if (isneg)
0095    {
0096       t = arg;
0097       t.negate();
0098       eval_exp(res, t);
0099       t.swap(res);
0100       res = limb_type(1);
0101       eval_divide(res, t);
0102       return;
0103    }
0104
0105    eval_divide(n, arg, default_ops::get_constant_ln2<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >());
0106    eval_floor(n, n);
0107    eval_multiply(t, n, default_ops::get_constant_ln2<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >());
0108    eval_subtract(t, arg);
0109    t.negate();
0110    if (t.compare(default_ops::get_constant_ln2<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >()) > 0)
0111    {
0112       // There are some rare cases where the multiply rounds down leaving a remainder > ln2
0113       // See https://github.com/boostorg/multiprecision/issues/120
0114       eval_increment(n);
0115       t = limb_type(0);
0116    }
0117    if (eval_get_sign(t) < 0)
0118    {
0119       // There are some very rare cases where arg/ln2 is an integer, and the subsequent multiply
0120       // rounds up, in that situation t ends up negative at this point which breaks our invariants below:
0121       t = limb_type(0);
0122    }
0123
0124    Exponent k, nn;
0125    eval_convert_to(&nn, n);
0126
0127    if (nn == (std::numeric_limits<Exponent>::max)())
0128    {
0129       // The result will necessarily oveflow:
0130       res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
0131       return;
0132    }
0133
0134    BOOST_MP_ASSERT(t.compare(default_ops::get_constant_ln2<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >()) < 0);
0135
0136    k = nn ? Exponent(1) << (msb(nn) / 2) : 0;
0137    k = (std::min)(k, (Exponent)(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count / 4));
0138    eval_ldexp(t, t, -k);
0139
0140    eval_exp_taylor(res, t);
0141    //
0142    // Square 1 + res k times:
0143    //
0144    for (Exponent s = 0; s < k; ++s)
0145    {
0146       t.swap(res);
0147       eval_multiply(res, t, t);
0148       eval_ldexp(t, t, 1);
0149       eval_add(res, t);
0150    }
0151    eval_add(res, limb_type(1));
0152    eval_ldexp(res, res, nn);
0153 }
0154
0155 }}} // namespace boost::multiprecision::backends
0156
0157 #endif