ceph/src/boost/boost/multiprecision/cpp_bin_float/transcendental.hpp

   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_MULTIPRECISION_CPP_BIN_FLOAT_TRANSCENDENTAL_HPP
   7 #define BOOST_MULTIPRECISION_CPP_BIN_FLOAT_TRANSCENDENTAL_HPP
   8
   9 #include <boost/multiprecision/detail/assert.hpp>
  10
  11 namespace boost { namespace multiprecision { namespace backends {
  12
  13 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  14 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)
  15 {
  16    constexpr const 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
  36 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  37 void 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