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