From 364a26b757f6a27239b2f9c9ac876c1e7a93f839 Mon Sep 17 00:00:00 2001
From: lyk <lyk-boya@outlook.com>
Date: Sat, 22 Jun 2024 20:13:15 +0800
Subject: [PATCH 1/7] Chebyshev approximation of sin(x) with degree 29

---
 include/gcem_incl/sin.hpp | 26 +++++++++++++++-----------
 tests/sin.cpp             |  5 ++++-
 2 files changed, 19 insertions(+), 12 deletions(-)

diff --git a/include/gcem_incl/sin.hpp b/include/gcem_incl/sin.hpp
index 8aa352d..3f224ca 100644
--- a/include/gcem_incl/sin.hpp
+++ b/include/gcem_incl/sin.hpp
@@ -39,6 +39,20 @@ noexcept
     return T(2)*x/(T(1) + x*x);
 }
 
+template<typename T>
+constexpr
+T
+sin_compute_chebyshev(const T x)
+noexcept {
+    constexpr T two_pie = 2 * T(GCEM_PI);
+    // limit x to range(-pi, pi)
+    const T x1 = x - two_pie * floor(x / two_pie) - T(GCEM_PI);
+    const T x2 = x1 * x1;
+    // Chebyshev approximation of sin(x) with degree 29
+    // weights are from: https://publik-void.github.io/sin-cos-approximations/#_sin_rel_error_minimized_degree_29
+    return -x1*(1. + x2*(-0.166666666666666666666666666666666629 + x2*(0.00833333333333333333333333333333222318 + x2*(-0.000198412698412698412698412698399723487 + x2*(2.75573192239858906525573184281026561e-6 + x2*(-2.50521083854417187750518138734779333e-8 + x2*(1.60590438368216145993211483359730933e-10 + x2*(-7.64716373181981646407639862566185217e-13 + x2*(2.81145725434551937513944561966331531e-15 + x2*(-8.22063524662315930664562316486611394e-18 + x2*(1.9572941062679701610307767721364446e-20 + x2*(-3.8681701397118386284141935821892017e-23 + x2*(6.44694091890808795934274904084719604e-26 + x2*(-9.18181103282828321298337267365314117e-29 + 1.10854165613066936680763885008011392e-31*x2))))))))))))));
+}
+
 template<typename T>
 constexpr
 T
@@ -51,17 +65,7 @@ noexcept
             // indistinguishable from zero
             GCLIM<T>::min() > abs(x) ? \
                 T(0) :
-            // special cases: pi/2 and pi
-            GCLIM<T>::min() > abs(x - T(GCEM_HALF_PI)) ? \
-                T(1) :
-            GCLIM<T>::min() > abs(x + T(GCEM_HALF_PI)) ? \
-                - T(1) :
-            GCLIM<T>::min() > abs(x - T(GCEM_PI)) ? \
-                T(0) :
-            GCLIM<T>::min() > abs(x + T(GCEM_PI)) ? \
-                - T(0) :
-            // else
-                sin_compute( tan(x/T(2)) ) );
+                sin_compute_chebyshev(x));
 }
 
 }
diff --git a/tests/sin.cpp b/tests/sin.cpp
index bf080f7..0f6f746 100644
--- a/tests/sin.cpp
+++ b/tests/sin.cpp
@@ -37,7 +37,10 @@ int main()
     GCEM_TEST_COMPARE_VALS(gcem::sin,std::sin,11.1L);
     GCEM_TEST_COMPARE_VALS(gcem::sin,std::sin,50.0L);
     GCEM_TEST_COMPARE_VALS(gcem::sin,std::sin,150.0L);
-
+    GCEM_TEST_COMPARE_VALS(gcem::sin,std::sin,GCEM_PI);
+    GCEM_TEST_COMPARE_VALS(gcem::sin,std::sin,GCEM_HALF_PI);
+    GCEM_TEST_COMPARE_VALS(gcem::sin,std::sin,-GCEM_PI);
+    GCEM_TEST_COMPARE_VALS(gcem::sin,std::sin,-GCEM_HALF_PI);
     GCEM_TEST_COMPARE_VALS(gcem::sin,std::sin,TEST_NAN);
 
     //

From 7c973d6f130f10c45ed5301547a4d0ecad11d411 Mon Sep 17 00:00:00 2001
From: lyk <lyk-boya@outlook.com>
Date: Sat, 22 Jun 2024 20:20:30 +0800
Subject: [PATCH 2/7] Chebyshev approximation of cos(x) with degree 30

---
 include/gcem_incl/cos.hpp | 26 +++++++++++++++-----------
 tests/cos.cpp             |  5 ++++-
 2 files changed, 19 insertions(+), 12 deletions(-)

diff --git a/include/gcem_incl/cos.hpp b/include/gcem_incl/cos.hpp
index 3f67460..7197fd5 100644
--- a/include/gcem_incl/cos.hpp
+++ b/include/gcem_incl/cos.hpp
@@ -37,6 +37,20 @@ noexcept
     return( T(1) - x*x)/(T(1) + x*x );
 }
 
+template<typename T>
+constexpr
+T
+cos_compute_chebyshev(const T x)
+noexcept {
+    constexpr T two_pie = 2 * T(GCEM_PI);
+    // limit x to range(0, 2*pi)
+    const T x1 = x - two_pie * floor(x / two_pie);
+    const T x2 = x1 * x1;
+    // Chebyshev approximation of sin(x) with degree 30
+    // weights are from: https://publik-void.github.io/sin-cos-approximations/#_cos_rel_error_minimized_degree_30
+    return 1. + x2*(-0.499999999999999999999999999999999997 + x2*(0.0416666666666666666666666666666665864 + x2*(-0.00138888888888888888888888888888792216 + x2*(0.0000248015873015873015873015872954247692 + x2*(-2.75573192239858906525573168323374212e-7 + x2*(2.08767569878680989792094774070130099e-9 + x2*(-1.147074559772972471374259582051505e-11 + x2*(4.77947733238738528351115798998742201e-14 + x2*(-1.56192069685862134697357839719554288e-16 + x2*(4.11031762331127151219297331433042247e-19 + x2*(-8.8967913919984472368244650875181655e-22 + x2*(1.61173755446230144849643763129805917e-24 + x2*(-2.47959193649582299714152529144837331e-27 + x2*(3.27913485904815466449434894827885198e-30 - 3.69081529290172836989264646507263976e-33*x2))))))))))))));
+}
+
 template<typename T>
 constexpr
 T
@@ -49,17 +63,7 @@ noexcept
             // indistinguishable from 0
             GCLIM<T>::min() > abs(x) ? 
                 T(1) :
-            // special cases: pi/2 and pi
-            GCLIM<T>::min() > abs(x - T(GCEM_HALF_PI)) ? \
-                T(0) :
-            GCLIM<T>::min() > abs(x + T(GCEM_HALF_PI)) ? \
-                T(0) :
-            GCLIM<T>::min() > abs(x - T(GCEM_PI)) ? \
-                - T(1) :
-            GCLIM<T>::min() > abs(x + T(GCEM_PI)) ? \
-                - T(1) :
-            // else
-                cos_compute( tan(x/T(2)) ) );
+                cos_compute_chebyshev(x));
 }
 
 }
diff --git a/tests/cos.cpp b/tests/cos.cpp
index 27ac1be..1ff4c6b 100644
--- a/tests/cos.cpp
+++ b/tests/cos.cpp
@@ -37,7 +37,10 @@ int main()
     GCEM_TEST_COMPARE_VALS(gcem::cos,std::cos,11.1L);
     GCEM_TEST_COMPARE_VALS(gcem::cos,std::cos,50.0L);
     GCEM_TEST_COMPARE_VALS(gcem::cos,std::cos,150.0L);
-
+    GCEM_TEST_COMPARE_VALS(gcem::cos,std::cos,GCEM_PI);
+    GCEM_TEST_COMPARE_VALS(gcem::cos,std::cos,GCEM_HALF_PI);
+    GCEM_TEST_COMPARE_VALS(gcem::cos,std::cos,-GCEM_PI);
+    GCEM_TEST_COMPARE_VALS(gcem::cos,std::cos,-GCEM_HALF_PI);
     GCEM_TEST_COMPARE_VALS(gcem::cos,std::cos,TEST_NAN);
 
     //

From f7de7ce07833c2cfcbedc767a5e9fa14f716c5cb Mon Sep 17 00:00:00 2001
From: lyk <lyk-boya@outlook.com>
Date: Sat, 22 Jun 2024 20:37:31 +0800
Subject: [PATCH 3/7] optimize tan using fast sin and fast cos

---
 include/gcem.hpp          |  3 ++-
 include/gcem_incl/cos.hpp |  2 +-
 include/gcem_incl/sin.hpp |  4 +---
 include/gcem_incl/tan.hpp | 20 +++++++++++++++++---
 4 files changed, 21 insertions(+), 8 deletions(-)

diff --git a/include/gcem.hpp b/include/gcem.hpp
index 4da27db..f66eb1e 100644
--- a/include/gcem.hpp
+++ b/include/gcem.hpp
@@ -70,9 +70,10 @@ namespace gcem
     #include "gcem_incl/gcd.hpp"
     #include "gcem_incl/lcm.hpp"
 
-    #include "gcem_incl/tan.hpp"
     #include "gcem_incl/cos.hpp"
     #include "gcem_incl/sin.hpp"
+    #include "gcem_incl/tan.hpp"
+
 
     #include "gcem_incl/atan.hpp"
     #include "gcem_incl/atan2.hpp"
diff --git a/include/gcem_incl/cos.hpp b/include/gcem_incl/cos.hpp
index 7197fd5..f5b1693 100644
--- a/include/gcem_incl/cos.hpp
+++ b/include/gcem_incl/cos.hpp
@@ -19,7 +19,7 @@
   ################################################################################*/
 
 /*
- * compile-time cosine function using tan(x/2)
+ * compile-time cosine function using 30th degree Chebyshev approximation
  */
 
 #ifndef _gcem_cos_HPP
diff --git a/include/gcem_incl/sin.hpp b/include/gcem_incl/sin.hpp
index 3f224ca..764a11e 100644
--- a/include/gcem_incl/sin.hpp
+++ b/include/gcem_incl/sin.hpp
@@ -19,9 +19,7 @@
   ################################################################################*/
 
 /*
- * compile-time sine function using tan(x/2)
- * 
- * see eq. 5.4.8 in Numerical Recipes
+ * compile-time sine function using 29th degree Chebyshev polynomial approximation
  */
 
 #ifndef _gcem_sin_HPP
diff --git a/include/gcem_incl/tan.hpp b/include/gcem_incl/tan.hpp
index 8872570..8d95fb2 100644
--- a/include/gcem_incl/tan.hpp
+++ b/include/gcem_incl/tan.hpp
@@ -97,6 +97,19 @@ noexcept
                 tan_cf_main(x) );
 }
 
+template<typename T>
+constexpr
+T
+tan_from_sin_cos(const T x)
+noexcept
+{
+    // There should be faster ways to calculate tan(x) directly
+    // At least as fast as calculating sin(x) and cos(x) using Chebyshev polynomials
+    const T sin_x = sin(x);
+    const T cos_x = cos(x);
+    return( sin_x/cos_x );
+}
+
 template<typename T>
 constexpr
 T
@@ -110,9 +123,10 @@ noexcept
             GCLIM<T>::min() > abs(x) ? \
                 T(0) :
             // else
-                x < T(0) ? \
-                    - tan_begin(-x) : 
-                      tan_begin( x) );
+                tan_from_sin_cos(x));
+                // x < T(0) ? \
+                //     - tan_begin(-x) :
+                //       tan_begin( x) );
 }
 
 }

From d7897c99bd583fc6c29444b764005b54edb5af9e Mon Sep 17 00:00:00 2001
From: lyk <lyk-boya@outlook.com>
Date: Sat, 22 Jun 2024 21:17:02 +0800
Subject: [PATCH 4/7] deal with tan special cases (pi/2)

---
 include/gcem_incl/tan.hpp | 10 +++++++---
 tests/tan.cpp             |  6 +++++-
 2 files changed, 12 insertions(+), 4 deletions(-)

diff --git a/include/gcem_incl/tan.hpp b/include/gcem_incl/tan.hpp
index 8d95fb2..629bc31 100644
--- a/include/gcem_incl/tan.hpp
+++ b/include/gcem_incl/tan.hpp
@@ -105,9 +105,13 @@ noexcept
 {
     // There should be faster ways to calculate tan(x) directly
     // At least as fast as calculating sin(x) and cos(x) using Chebyshev polynomials
-    const T sin_x = sin(x);
-    const T cos_x = cos(x);
-    return( sin_x/cos_x );
+    // limit x to range(0, pi)
+    T x1 = x - T(GCEM_PI) * floor(x / T(GCEM_PI));
+    // calculate
+    if(x1 < T(GCEM_HALF_PI))    return sin(x1) / cos(x1);
+    if(x1 == T(GCEM_HALF_PI))   return GCLIM<T>::quiet_NaN();
+    x1 = T(GCEM_PI) - x1;
+    return - sin(x1) / cos(x1);
 }
 
 template<typename T>
diff --git a/tests/tan.cpp b/tests/tan.cpp
index 291cb65..75c0770 100644
--- a/tests/tan.cpp
+++ b/tests/tan.cpp
@@ -31,6 +31,7 @@ int main()
 
     // static constexpr long double lrgval = std::numeric_limits<int>::max()*100.0L;
 
+    std::cout << "tan(pi/2) = " << std::tan(GCEM_HALF_PI) << std::endl;
     GCEM_TEST_COMPARE_VALS(gcem::tan,std::tan, 0.0L);
     GCEM_TEST_COMPARE_VALS(gcem::tan,std::tan, 0.001L);
     GCEM_TEST_COMPARE_VALS(gcem::tan,std::tan, 1.001L);
@@ -39,7 +40,10 @@ int main()
     GCEM_TEST_COMPARE_VALS(gcem::tan,std::tan, 50.0L);
     GCEM_TEST_COMPARE_VALS(gcem::tan,std::tan, -1.5L);
     // GCEM_TEST_COMPARE_VALS(gcem::tan,std::tan, lrgval);
-
+    GCEM_TEST_COMPARE_VALS(gcem::tan,std::tan, GCEM_HALF_PI - 1e-6);
+    GCEM_TEST_COMPARE_VALS(gcem::tan,std::tan, GCEM_HALF_PI + 1e-6);
+    GCEM_TEST_COMPARE_VALS(gcem::tan,std::tan, -GCEM_HALF_PI - 1e-6);
+    GCEM_TEST_COMPARE_VALS(gcem::tan,std::tan, -GCEM_HALF_PI + 1e-6);
     GCEM_TEST_COMPARE_VALS(gcem::tan,std::tan, TEST_NAN);
 
     //

From 2e58648e832b79bc1f2a3f6694ecb075adbd36cf Mon Sep 17 00:00:00 2001
From: lyk <lyk-boya@outlook.com>
Date: Sat, 22 Jun 2024 21:54:21 +0800
Subject: [PATCH 5/7] optimize sqrt using Babylonian Method

---
 include/gcem_incl/sqrt.hpp | 18 +++++++++++++++++-
 1 file changed, 17 insertions(+), 1 deletion(-)

diff --git a/include/gcem_incl/sqrt.hpp b/include/gcem_incl/sqrt.hpp
index 4d5555d..20abdc2 100644
--- a/include/gcem_incl/sqrt.hpp
+++ b/include/gcem_incl/sqrt.hpp
@@ -64,6 +64,21 @@ noexcept
                 m_val * sqrt_recur(x, x / T(2), 0) );
 }
 
+template<typename T>
+T
+sqrt_babylonian(const T x)
+noexcept
+{
+    // see here for more details: https://blogs.sas.com/content/iml/2016/05/16/babylonian-square-roots.html
+    T xn = x;
+    T xn1 = (xn + x/xn) / T(2);
+    while(abs(xn1 - xn) > GCLIM<T>::epsilon()) {
+        xn = xn1;
+        xn1 = (xn + x/xn) / T(2);
+    }
+    return xn1;
+}
+
 template<typename T>
 constexpr
 T
@@ -84,7 +99,8 @@ noexcept
             GCLIM<T>::min() > abs(T(1) - x) ? \
                 x :
             // else
-            sqrt_simplify(x, T(1)) );
+            // sqrt_simplify(x, T(1)) );
+            sqrt_babylonian(x) );
 }
 
 }

From e338f6c0074eaf86e670853dfeed0a081dbac791 Mon Sep 17 00:00:00 2001
From: lyk <lyk-boya@outlook.com>
Date: Sat, 22 Jun 2024 22:04:41 +0800
Subject: [PATCH 6/7] add GCEM_2PI

---
 include/gcem_incl/cos.hpp          | 3 +--
 include/gcem_incl/gcem_options.hpp | 4 ++++
 include/gcem_incl/sin.hpp          | 3 +--
 3 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/include/gcem_incl/cos.hpp b/include/gcem_incl/cos.hpp
index f5b1693..99da2cc 100644
--- a/include/gcem_incl/cos.hpp
+++ b/include/gcem_incl/cos.hpp
@@ -42,9 +42,8 @@ constexpr
 T
 cos_compute_chebyshev(const T x)
 noexcept {
-    constexpr T two_pie = 2 * T(GCEM_PI);
     // limit x to range(0, 2*pi)
-    const T x1 = x - two_pie * floor(x / two_pie);
+    const T x1 = x - T(GCEM_2PI) * floor(x / T(GCEM_2PI));
     const T x2 = x1 * x1;
     // Chebyshev approximation of sin(x) with degree 30
     // weights are from: https://publik-void.github.io/sin-cos-approximations/#_cos_rel_error_minimized_degree_30
diff --git a/include/gcem_incl/gcem_options.hpp b/include/gcem_incl/gcem_options.hpp
index 1e988c3..2762abf 100644
--- a/include/gcem_incl/gcem_options.hpp
+++ b/include/gcem_incl/gcem_options.hpp
@@ -98,6 +98,10 @@ namespace gcem
     #define GCEM_PI 3.1415926535897932384626433832795028841972L
 #endif
 
+#ifndef GCEM_2PI
+    #define GCEM_2PI 6.2831853071795864769252867665590057683943L
+#endif
+
 #ifndef GCEM_LOG_PI
     #define GCEM_LOG_PI 1.1447298858494001741434273513530587116473L
 #endif
diff --git a/include/gcem_incl/sin.hpp b/include/gcem_incl/sin.hpp
index 764a11e..e48bc35 100644
--- a/include/gcem_incl/sin.hpp
+++ b/include/gcem_incl/sin.hpp
@@ -42,9 +42,8 @@ constexpr
 T
 sin_compute_chebyshev(const T x)
 noexcept {
-    constexpr T two_pie = 2 * T(GCEM_PI);
     // limit x to range(-pi, pi)
-    const T x1 = x - two_pie * floor(x / two_pie) - T(GCEM_PI);
+    const T x1 = x - T(GCEM_2PI) * floor(x / T(GCEM_2PI)) - T(GCEM_PI);
     const T x2 = x1 * x1;
     // Chebyshev approximation of sin(x) with degree 29
     // weights are from: https://publik-void.github.io/sin-cos-approximations/#_sin_rel_error_minimized_degree_29

From 7d8f81390e6e8b8ae6f45ab711131a8c92fa6c84 Mon Sep 17 00:00:00 2001
From: lyk <lyk-boya@outlook.com>
Date: Sat, 22 Jun 2024 22:53:17 +0800
Subject: [PATCH 7/7] remove declaration in 'constexpr' function body

---
 include/gcem_incl/cos.hpp  | 15 ++++++++++-----
 include/gcem_incl/sin.hpp  | 15 ++++++++++-----
 include/gcem_incl/sqrt.hpp | 30 ++++++++++++++++++++----------
 include/gcem_incl/tan.hpp  | 23 ++++++++++++++++-------
 4 files changed, 56 insertions(+), 27 deletions(-)

diff --git a/include/gcem_incl/cos.hpp b/include/gcem_incl/cos.hpp
index 99da2cc..de8187e 100644
--- a/include/gcem_incl/cos.hpp
+++ b/include/gcem_incl/cos.hpp
@@ -40,16 +40,21 @@ noexcept
 template<typename T>
 constexpr
 T
-cos_compute_chebyshev(const T x)
+cos_compute_chebyshev_inner(const T x2)
 noexcept {
-    // limit x to range(0, 2*pi)
-    const T x1 = x - T(GCEM_2PI) * floor(x / T(GCEM_2PI));
-    const T x2 = x1 * x1;
     // Chebyshev approximation of sin(x) with degree 30
     // weights are from: https://publik-void.github.io/sin-cos-approximations/#_cos_rel_error_minimized_degree_30
     return 1. + x2*(-0.499999999999999999999999999999999997 + x2*(0.0416666666666666666666666666666665864 + x2*(-0.00138888888888888888888888888888792216 + x2*(0.0000248015873015873015873015872954247692 + x2*(-2.75573192239858906525573168323374212e-7 + x2*(2.08767569878680989792094774070130099e-9 + x2*(-1.147074559772972471374259582051505e-11 + x2*(4.77947733238738528351115798998742201e-14 + x2*(-1.56192069685862134697357839719554288e-16 + x2*(4.11031762331127151219297331433042247e-19 + x2*(-8.8967913919984472368244650875181655e-22 + x2*(1.61173755446230144849643763129805917e-24 + x2*(-2.47959193649582299714152529144837331e-27 + x2*(3.27913485904815466449434894827885198e-30 - 3.69081529290172836989264646507263976e-33*x2))))))))))))));
 }
 
+template<typename T>
+constexpr
+T
+cos_compute_chebyshev(const T x1)
+noexcept {
+    return cos_compute_chebyshev_inner(x1 * x1);
+}
+
 template<typename T>
 constexpr
 T
@@ -62,7 +67,7 @@ noexcept
             // indistinguishable from 0
             GCLIM<T>::min() > abs(x) ? 
                 T(1) :
-                cos_compute_chebyshev(x));
+                cos_compute_chebyshev(x - T(GCEM_2PI) * floor(x / T(GCEM_2PI))));
 }
 
 }
diff --git a/include/gcem_incl/sin.hpp b/include/gcem_incl/sin.hpp
index e48bc35..f980cde 100644
--- a/include/gcem_incl/sin.hpp
+++ b/include/gcem_incl/sin.hpp
@@ -40,16 +40,21 @@ noexcept
 template<typename T>
 constexpr
 T
-sin_compute_chebyshev(const T x)
+sin_compute_chebyshev_inner(const T x1, const T x2)
 noexcept {
-    // limit x to range(-pi, pi)
-    const T x1 = x - T(GCEM_2PI) * floor(x / T(GCEM_2PI)) - T(GCEM_PI);
-    const T x2 = x1 * x1;
     // Chebyshev approximation of sin(x) with degree 29
     // weights are from: https://publik-void.github.io/sin-cos-approximations/#_sin_rel_error_minimized_degree_29
     return -x1*(1. + x2*(-0.166666666666666666666666666666666629 + x2*(0.00833333333333333333333333333333222318 + x2*(-0.000198412698412698412698412698399723487 + x2*(2.75573192239858906525573184281026561e-6 + x2*(-2.50521083854417187750518138734779333e-8 + x2*(1.60590438368216145993211483359730933e-10 + x2*(-7.64716373181981646407639862566185217e-13 + x2*(2.81145725434551937513944561966331531e-15 + x2*(-8.22063524662315930664562316486611394e-18 + x2*(1.9572941062679701610307767721364446e-20 + x2*(-3.8681701397118386284141935821892017e-23 + x2*(6.44694091890808795934274904084719604e-26 + x2*(-9.18181103282828321298337267365314117e-29 + 1.10854165613066936680763885008011392e-31*x2))))))))))))));
 }
 
+template<typename T>
+constexpr
+T
+sin_compute_chebyshev(const T x1)
+noexcept {
+    return sin_compute_chebyshev_inner(x1, x1 * x1);
+}
+
 template<typename T>
 constexpr
 T
@@ -62,7 +67,7 @@ noexcept
             // indistinguishable from zero
             GCLIM<T>::min() > abs(x) ? \
                 T(0) :
-                sin_compute_chebyshev(x));
+                sin_compute_chebyshev(x - T(GCEM_2PI) * floor(x / T(GCEM_2PI)) - T(GCEM_PI)));
 }
 
 }
diff --git a/include/gcem_incl/sqrt.hpp b/include/gcem_incl/sqrt.hpp
index 20abdc2..d171ce1 100644
--- a/include/gcem_incl/sqrt.hpp
+++ b/include/gcem_incl/sqrt.hpp
@@ -64,19 +64,29 @@ noexcept
                 m_val * sqrt_recur(x, x / T(2), 0) );
 }
 
+// template<typename T>
+// T
+// sqrt_babylonian(const T x)
+// noexcept
+// {
+//     // see here for more details: https://blogs.sas.com/content/iml/2016/05/16/babylonian-square-roots.html
+//     T xn = x;
+//     T xn1 = (xn + x/xn) / T(2);
+//     while(abs(xn1 - xn) > GCLIM<T>::epsilon()) {
+//         xn = xn1;
+//         xn1 = (xn + x/xn) / T(2);
+//     }
+//     return xn1;
+// }
+
 template<typename T>
 T
-sqrt_babylonian(const T x)
+sqrt_babylonian(const T x, const T xn, const T xn1)
 noexcept
 {
-    // see here for more details: https://blogs.sas.com/content/iml/2016/05/16/babylonian-square-roots.html
-    T xn = x;
-    T xn1 = (xn + x/xn) / T(2);
-    while(abs(xn1 - xn) > GCLIM<T>::epsilon()) {
-        xn = xn1;
-        xn1 = (xn + x/xn) / T(2);
-    }
-    return xn1;
+    return ( abs(xn1 - xn) < GCLIM<T>::epsilon() ? \
+                xn1 :
+                sqrt_babylonian(x, xn1, (xn1 + x/xn1) / T(2)));
 }
 
 template<typename T>
@@ -100,7 +110,7 @@ noexcept
                 x :
             // else
             // sqrt_simplify(x, T(1)) );
-            sqrt_babylonian(x) );
+            sqrt_babylonian(x, x, (x + T(1)) / T(2)) );
 }
 
 }
diff --git a/include/gcem_incl/tan.hpp b/include/gcem_incl/tan.hpp
index 629bc31..ab00b8b 100644
--- a/include/gcem_incl/tan.hpp
+++ b/include/gcem_incl/tan.hpp
@@ -100,18 +100,27 @@ noexcept
 template<typename T>
 constexpr
 T
-tan_from_sin_cos(const T x)
+tan_from_sin_cos_inner(const T x1)
 noexcept
 {
     // There should be faster ways to calculate tan(x) directly
     // At least as fast as calculating sin(x) and cos(x) using Chebyshev polynomials
+    return (
+        x1 < T(GCEM_HALF_PI) ? sin(x1) / cos(x1) :
+        x1 > T(GCEM_HALF_PI) ? - sin(T(GCEM_PI) - x1) / cos(T(GCEM_PI) - x1) :
+        GCLIM<T>::quiet_NaN()   // x1 == 1/2 pi
+    );
+}
+
+
+template<typename T>
+constexpr
+T
+tan_from_sin_cos(const T x)
+noexcept
+{
     // limit x to range(0, pi)
-    T x1 = x - T(GCEM_PI) * floor(x / T(GCEM_PI));
-    // calculate
-    if(x1 < T(GCEM_HALF_PI))    return sin(x1) / cos(x1);
-    if(x1 == T(GCEM_HALF_PI))   return GCLIM<T>::quiet_NaN();
-    x1 = T(GCEM_PI) - x1;
-    return - sin(x1) / cos(x1);
+    return tan_from_sin_cos_inner(x - T(GCEM_PI) * floor(x / T(GCEM_PI)));
 }
 
 template<typename T>