From 8422f5307b0498d09cf626f38acb03fe9f1efd94 Mon Sep 17 00:00:00 2001
From: Keith O'Hara <kth.ohr@gmail.com>
Date: Thu, 2 Mar 2023 17:49:12 -0500
Subject: [PATCH] add C++14 implementation of exp

---
 include/gcem_incl/exp.hpp | 34 +++++++++++++++++++++++++++++-----
 1 file changed, 29 insertions(+), 5 deletions(-)

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<typename T>
+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<typename T>
 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<typename T>
 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<typename T>
@@ -72,7 +96,7 @@ noexcept
             //
             is_neginf(x) ? \
                 T(0) :
-            //
+            // indistinguishable from zero
             GCLIM<T>::min() > abs(x) ? \
                 T(1) : 
             //