diff --git a/include/gcem_incl/exp.hpp b/include/gcem_incl/exp.hpp index 7fd63fe..595ffc2 100644 --- a/include/gcem_incl/exp.hpp +++ b/include/gcem_incl/exp.hpp @@ -28,6 +28,30 @@ namespace internal { +// see https://en.wikipedia.org/wiki/Euler%27s_continued_fraction_formula + +#if __cplusplus >= 201402L // C++14 version + +template +constexpr +T +exp_cf_recur(const T x, const int depth_end) +noexcept +{ + int depth = GCEM_EXP_MAX_ITER_SMALL - 1; + T res = T(1); + + while (depth > depth_end - 1) { + res = T(1) + x/T(depth - 1) - x/depth/res; + + --depth; + } + + return res; +} + +#else // C++11 version + template constexpr T @@ -36,20 +60,20 @@ noexcept { return( depth < GCEM_EXP_MAX_ITER_SMALL ? \ // if - depth == 1 ? \ - T(1) - x/exp_cf_recur(x,depth+1) : - T(1) + x/T(depth - 1) - x/depth/exp_cf_recur(x,depth+1) : + T(1) + x/T(depth - 1) - x/depth/exp_cf_recur(x,depth+1) : // else T(1) ); } +#endif + template constexpr T exp_cf(const T x) noexcept { - return( T(1)/exp_cf_recur(x,1) ); + return( T(1) / (T(1) - x / exp_cf_recur(x,2)) ); } template @@ -72,7 +96,7 @@ noexcept // is_neginf(x) ? \ T(0) : - // + // indistinguishable from zero GCLIM::min() > abs(x) ? \ T(1) : //