Chris@102: ///////////////////////////////////////////////////////////////
Chris@102: //  Copyright 2013 John Maddock. Distributed under the Boost
Chris@102: //  Software License, Version 1.0. (See accompanying file
Chris@102: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_
Chris@102: 
Chris@102: #ifndef BOOST_MULTIPRECISION_CPP_BIN_FLOAT_TRANSCENDENTAL_HPP
Chris@102: #define BOOST_MULTIPRECISION_CPP_BIN_FLOAT_TRANSCENDENTAL_HPP
Chris@102: 
Chris@102: namespace boost{ namespace multiprecision{ namespace backends{
Chris@102: 
Chris@102: template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
Chris@102: 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)
Chris@102: {
Chris@102:    static const int bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count;
Chris@102:    //
Chris@102:    // Taylor series for small argument, note returns exp(x) - 1:
Chris@102:    //
Chris@102:    res = limb_type(0);
Chris@102:    cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> num(arg), denom, t;
Chris@102:    denom = limb_type(1);
Chris@102:    eval_add(res, num);
Chris@102: 
Chris@102:    for(unsigned k = 2; ; ++k)
Chris@102:    {
Chris@102:       eval_multiply(denom, k);
Chris@102:       eval_multiply(num, arg);
Chris@102:       eval_divide(t, num, denom);
Chris@102:       eval_add(res, t);
Chris@102:       if(eval_is_zero(t) || (res.exponent() - bits > t.exponent()))
Chris@102:          break;
Chris@102:    }
Chris@102: }
Chris@102: 
Chris@102: template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
Chris@102: void eval_exp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
Chris@102: {
Chris@102:    //
Chris@102:    // This is based on MPFR's method, let:
Chris@102:    //
Chris@102:    // n = floor(x / ln(2))
Chris@102:    //
Chris@102:    // Then:
Chris@102:    //
Chris@102:    // r = x - n ln(2) : 0 <= r < ln(2)
Chris@102:    //
Chris@102:    // We can reduce r further by dividing by 2^k, with k ~ sqrt(n),
Chris@102:    // so if:
Chris@102:    //
Chris@102:    // e0 = exp(r / 2^k) - 1
Chris@102:    //
Chris@102:    // With e0 evaluated by taylor series for small arguments, then:
Chris@102:    //
Chris@102:    // exp(x) = 2^n (1 + e0)^2^k
Chris@102:    //
Chris@102:    // Note that to preserve precision we actually square (1 + e0) k times, calculating
Chris@102:    // the result less one each time, i.e.
Chris@102:    //
Chris@102:    // (1 + e0)^2 - 1 = e0^2 + 2e0
Chris@102:    //
Chris@102:    // Then add the final 1 at the end, given that e0 is small, this effectively wipes
Chris@102:    // out the error in the last step.
Chris@102:    //
Chris@102:    using default_ops::eval_multiply;
Chris@102:    using default_ops::eval_subtract;
Chris@102:    using default_ops::eval_add;
Chris@102:    using default_ops::eval_convert_to;
Chris@102: 
Chris@102:    int type = eval_fpclassify(arg);
Chris@102:    bool isneg = eval_get_sign(arg) < 0;
Chris@102:    if(type == (int)FP_NAN)
Chris@102:    {
Chris@102:       res = arg;
Chris@102:       return;
Chris@102:    }
Chris@102:    else if(type == (int)FP_INFINITE)
Chris@102:    {
Chris@102:       res = arg;
Chris@102:       if(isneg)
Chris@102:          res = limb_type(0u);
Chris@102:       else 
Chris@102:          res = arg;
Chris@102:       return;
Chris@102:    }
Chris@102:    else if(type == (int)FP_ZERO)
Chris@102:    {
Chris@102:       res = limb_type(1);
Chris@102:       return;
Chris@102:    }
Chris@102:    cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> t, n;
Chris@102:    if(isneg)
Chris@102:    {
Chris@102:       t = arg;
Chris@102:       t.negate();
Chris@102:       eval_exp(res, t);
Chris@102:       t.swap(res);
Chris@102:       res = limb_type(1);
Chris@102:       eval_divide(res, t);
Chris@102:       return;
Chris@102:    }
Chris@102: 
Chris@102:    eval_divide(n, arg, default_ops::get_constant_ln2<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >());
Chris@102:    eval_floor(n, n);
Chris@102:    eval_multiply(t, n, default_ops::get_constant_ln2<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >());
Chris@102:    eval_subtract(t, arg);
Chris@102:    t.negate();
Chris@102:    if(eval_get_sign(t) < 0)
Chris@102:    {
Chris@102:       // There are some very rare cases where arg/ln2 is an integer, and the subsequent multiply
Chris@102:       // rounds up, in that situation t ends up negative at this point which breaks our invariants below:
Chris@102:       t = limb_type(0);
Chris@102:    }
Chris@102:    BOOST_ASSERT(t.compare(default_ops::get_constant_ln2<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >()) < 0);
Chris@102: 
Chris@102:    Exponent k, nn;
Chris@102:    eval_convert_to(&nn, n);
Chris@102:    k = nn ? Exponent(1) << (msb(nn) / 2) : 0;
Chris@102:    eval_ldexp(t, t, -k);
Chris@102: 
Chris@102:    eval_exp_taylor(res, t);
Chris@102:    //
Chris@102:    // Square 1 + res k times:
Chris@102:    //
Chris@102:    for(int s = 0; s < k; ++s)
Chris@102:    {
Chris@102:       t.swap(res);
Chris@102:       eval_multiply(res, t, t);
Chris@102:       eval_ldexp(t, t, 1);
Chris@102:       eval_add(res, t);
Chris@102:    }
Chris@102:    eval_add(res, limb_type(1));
Chris@102:    eval_ldexp(res, res, nn);
Chris@102: }
Chris@102: 
Chris@102: }}} // namespaces
Chris@102: 
Chris@102: #endif
Chris@102: