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