diff --git a/doc/markdown/gamma.md b/doc/markdown/gamma.md
new file mode 100644
index 0000000..dcc3039
--- /dev/null
+++ b/doc/markdown/gamma.md
@@ -0,0 +1,18 @@
+# Gamma Functions
+
+## Airy Functions
+
+* **fac**,    [Factorial function](doubldoc.md#fac)
+* **gamma**,  [Gamma function](doubldoc.md#gamma)
+* **lgam**,   [Natural logarithm of gamma function](doubldoc.md#lgam)
+* **igam**,   [Incomplete gamma integral](doubldoc.md#igam)
+* **igamc**,  [Complemented incomplete gamma integral](doubldoc.md#igamc)
+* **igami**,  [Inverse of complemented imcomplete gamma integral](doubldoc.md#igami)
+* **incbet**, [Incomplete beta integral](doubldoc.md#incbet)
+* **incbi**,  [Inverse of imcomplete beta integral](doubldoc.md#incbi)
+* **psi**,    [Psi (digamma) function](doubldoc.md#psi)
+* **rgamma**, [Reciprocal gamma function](doubldoc.md#rgamma)
+
+## Beta Functions
+
+* **beta**,   [Beta function](doubldoc.md#beta)
diff --git a/include/cephes/bessel.h b/include/cephes/bessel.h
index 72a0425..831faec 100644
--- a/include/cephes/bessel.h
+++ b/include/cephes/bessel.h
@@ -5,6 +5,7 @@
  */
 #define CEPHES_BESSEL_H
 
+namespace cephes {
 #if defined(__cplusplus)
 extern "C" {
 #endif
@@ -46,6 +47,8 @@ double psi(double x);
 double struve(double v, double x);
 
 #if defined(__cplusplus)
-}
+} // extern "C"
 #endif
+}; // ::cephes
+
 #endif // CEPHES_BESSEL_H
\ No newline at end of file
diff --git a/include/cephes/exp_int.h b/include/cephes/exp_int.h
index 54b3009..3d1a6d8 100644
--- a/include/cephes/exp_int.h
+++ b/include/cephes/exp_int.h
@@ -5,6 +5,7 @@
  */
 #define CEPHES_EXP_INT_H
 
+namespace cephes {
 #if defined(__cplusplus)
 extern "C" {
 #endif
@@ -20,6 +21,7 @@ int sici(double x, double *si, double *ci);
 int shichi(double x, double *si, double *ci);
 
 #if defined(__cplusplus)
-}
+} // extern "C"
 #endif
+}; // ::cephes
 #endif // CEPHES_EXP_INT_H
\ No newline at end of file
diff --git a/include/cephes/gamma.h b/include/cephes/gamma.h
new file mode 100644
index 0000000..eab8851
--- /dev/null
+++ b/include/cephes/gamma.h
@@ -0,0 +1,44 @@
+#ifndef CEPHES_GAMMA_H
+/** Cephes double precision special functions suite
+ * 
+ * cephes/bessel
+ */
+#define CEPHES_GAMMA_H
+
+namespace cephes {
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+/** Gamma Functions */
+/* misc/fac.c */
+double fac(int i);
+
+/* cprob/gamma.c */
+double gamma(double x);
+double lgam(double x);
+/* cprob/igam.c */
+double igam(double a, double x);
+double igamc(double a, double x);
+/* cprob/igami.c */
+double igami(double a, double y0_);
+/* cprob/incbet.c */
+double incbet(double aa, double bb, double xx);
+/* cprob/incbi.c */
+double incbi(double aa, double bb, double yy0);
+
+/* misc/psi.c */
+double psi(double x);
+/* misc/rgamma.c */
+double rgamma(double x);
+
+/** Beta Functions */
+/* misc/beta.c */
+double beta(double a, double b);
+
+#if defined(__cplusplus)
+} // extern "C"
+#endif
+}; // ::cephes
+
+#endif // CEPHES_GAMMA_H
\ No newline at end of file
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index ec7a0e6..8b353d1 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -22,6 +22,6 @@ function(add_gtest test_name)
     gtest_discover_tests(${test_name})
 endfunction()
 
-add_subdirectory(bessel)
-add_gtest(airy)
 add_gtest(exp_int)
+add_subdirectory(gamma)
+add_subdirectory(bessel)
diff --git a/tests/bessel/CMakeLists.txt b/tests/bessel/CMakeLists.txt
index b56cf6c..527f7b6 100644
--- a/tests/bessel/CMakeLists.txt
+++ b/tests/bessel/CMakeLists.txt
@@ -1,3 +1,26 @@
 
+add_gtest(airy)
+
 add_gtest(j0)
 add_gtest(j1)
+add_gtest(jn)
+add_gtest(jv)
+
+add_gtest(i0)
+add_gtest(i0e)
+add_gtest(i1)
+add_gtest(i1e)
+add_gtest(iv)
+
+add_gtest(y0)
+add_gtest(y1)
+add_gtest(yn)
+
+add_gtest(k0)
+add_gtest(k0e)
+add_gtest(k1)
+add_gtest(k1e)
+add_gtest(kn)
+
+add_gtest(psi)
+add_gtest(struve)
diff --git a/tests/airy.cpp b/tests/bessel/airy.cpp
similarity index 86%
rename from tests/airy.cpp
rename to tests/bessel/airy.cpp
index 5b323e5..c0b3671 100644
--- a/tests/airy.cpp
+++ b/tests/bessel/airy.cpp
@@ -1,7 +1,6 @@
-#include <limits>
-#include <complex>
-#include <cephes/bessel.h>
 #include <xtest.hpp>
+#include <cephes/bessel.h>
+
 
 TEST(Airy, BasicAssertions) {
     const double nan64 = std::numeric_limits<double>::quiet_NaN();
@@ -9,7 +8,7 @@ TEST(Airy, BasicAssertions) {
     double x, ai, aip, bi, bip;
 
     x = 26.0;
-    ret = airy(x, &ai, &aip, &bi, &bip);
+    ret = cephes::airy(x, &ai, &aip, &bi, &bip);
     EXPECT_EQ(ret, -1);
     EXPECT_EQ(ai, 0.0);
     EXPECT_EQ(aip, 0.0);
@@ -24,7 +23,7 @@ TEST(Airy, Branches) {
 
     // x < -2.09
     x = -3.0;
-    ret = airy(x, &ai, &aip, &bi, &bip);
+    ret = cephes::airy(x, &ai, &aip, &bi, &bip);
     EXPECT_EQ(ret, 0);
     /*
         x = -3.0;
@@ -42,7 +41,7 @@ TEST(Airy, Branches) {
 
     // x >= 2.09 (cbrt(9))
     x = 5.0;
-    ret = airy(x, &ai, &aip, &bi, &bip);
+    ret = cephes::airy(x, &ai, &aip, &bi, &bip);
     EXPECT_EQ(ret, 0);
     /*
         x = 5.0;
@@ -60,7 +59,7 @@ TEST(Airy, Branches) {
 
     // x > 8.3203353
     x = 8.5;
-    ret = airy(x, &ai, &aip, &bi, &bip);
+    ret = cephes::airy(x, &ai, &aip, &bi, &bip);
     EXPECT_EQ(ret, 0);
     /*
         x = 8.5;
@@ -82,9 +81,9 @@ TEST(Hyp2f1, Errors) {
 
     // c is a negative integer
     // hypdiv:
-    EXPECT_GT(hyp2f1(2.0, 3.0, -4.0, 5.0), 1e308);
+    EXPECT_GT(cephes::hyp2f1(2.0, 3.0, -4.0, 5.0), 1e308);
     // ax > 1.0
-    EXPECT_GT(hyp2f1(2.0, 3.0, 4.0, 2.0), 1e308);
+    EXPECT_GT(cephes::hyp2f1(2.0, 3.0, 4.0, 2.0), 1e308);
 }
 
 TEST(Hyp2f1, Branches) {
diff --git a/tests/bessel/i0.cpp b/tests/bessel/i0.cpp
new file mode 100644
index 0000000..457d2fe
--- /dev/null
+++ b/tests/bessel/i0.cpp
@@ -0,0 +1,20 @@
+#include <xtest.hpp>
+#include <cephes/bessel.h>
+
+
+TEST(BesselI0, Branches) {
+    EXPECT_REL_NEAR_F64(cephes::i0(0.0), 1.0);
+
+    // x < 0
+    EXPECT_REL_NEAR_F64(cephes::i0(-0.0), 1.0);
+    EXPECT_REL_NEAR_F64(cephes::i0(-10.0), 2815.716628466253);
+
+    // x <= 8.0
+    EXPECT_REL_NEAR_F64(cephes::i0(1.0), 1.266065877752008);
+    EXPECT_REL_NEAR_F64(cephes::i0(8.0), 427.5641157218048);
+
+    // x > 8.0
+    EXPECT_REL_NEAR_F64(cephes::i0(9.0),  1093.588354511374);
+    EXPECT_REL_NEAR_F64(cephes::i0(10.0), 2815.716628466253);
+    EXPECT_REL_NEAR_F64(cephes::i0(100.0), 1.073751707131074e42);
+}
diff --git a/tests/bessel/i0e.cpp b/tests/bessel/i0e.cpp
new file mode 100644
index 0000000..f49456f
--- /dev/null
+++ b/tests/bessel/i0e.cpp
@@ -0,0 +1,23 @@
+#include <xtest.hpp>
+#include <cephes/bessel.h>
+
+/**
+    Table[NumberForm[Exp[-Abs[x]]*BesselI[nv, x], 16],
+        {x, {0.0, -10., 1.0, 8.0, 9.0, 10.0, 100.}}, {nv, {0}}]
+ */
+TEST(BesselI0Exp, Branches) {
+    EXPECT_REL_NEAR_F64(cephes::i0e(0.0), 1.0);
+
+    // x < 0
+    EXPECT_REL_NEAR_F64(cephes::i0e(-0.0), 1.0);
+    EXPECT_REL_NEAR_F64(cephes::i0e(-10.0), 0.1278333371634286);
+
+    // x <= 8.0
+    EXPECT_REL_NEAR_F64(cephes::i0e(1.0), 0.4657596075936404);
+    EXPECT_REL_NEAR_F64(cephes::i0e(8.0), 0.1434317818568503);
+
+    // x > 8.0
+    EXPECT_REL_NEAR_F64(cephes::i0e(9.0),  0.134959524581723);
+    EXPECT_REL_NEAR_F64(cephes::i0e(10.0), 0.1278333371634286);
+    EXPECT_REL_NEAR_F64(cephes::i0e(100.0), 0.03994437929909669);
+}
diff --git a/tests/bessel/i1.cpp b/tests/bessel/i1.cpp
new file mode 100644
index 0000000..b8ce566
--- /dev/null
+++ b/tests/bessel/i1.cpp
@@ -0,0 +1,23 @@
+#include <xtest.hpp>
+#include <cephes/bessel.h>
+
+/*
+    Table[NumberForm[BesselI[1, x], 16],
+        {x, {0.0, -10., 1.0, 8.0, 9.0, 10.0, 100.}}]
+*/
+TEST(BesselI1, Branches) {
+    EXPECT_REL_NEAR_F64(cephes::i1(0.0), 0.0);
+
+    // x < 0
+    EXPECT_REL_NEAR_F64(cephes::i1(-0.0), 0.0);
+    EXPECT_REL_NEAR_F64(cephes::i1(-10.0), -2670.988303701254);
+
+    // x <= 8.0
+    EXPECT_REL_NEAR_F64(cephes::i1(1.0), 0.5651591039924851);
+    EXPECT_REL_NEAR_F64(cephes::i1(8.0), 399.8731367825601);
+
+    // x > 8.0
+    EXPECT_REL_NEAR_F64(cephes::i1(9.0),  1030.914722516956);
+    EXPECT_REL_NEAR_F64(cephes::i1(10.0), 2670.988303701254);
+    EXPECT_REL_NEAR_F64(cephes::i1(100.0), 1.068369390338163e+42);
+}
diff --git a/tests/bessel/i1e.cpp b/tests/bessel/i1e.cpp
new file mode 100644
index 0000000..0a3ef46
--- /dev/null
+++ b/tests/bessel/i1e.cpp
@@ -0,0 +1,23 @@
+#include <xtest.hpp>
+#include <cephes/bessel.h>
+
+/**
+    Table[NumberForm[Exp[-Abs[x]]*BesselI[1, x], 16],
+        {x, {0.0, -10., 1.0, 8.0, 9.0, 10.0, 100.}}]
+ */
+TEST(BesselI1Exp, Branches) {
+    EXPECT_REL_NEAR_F64(cephes::i1e(0.0), 0.0);
+
+    // x < 0
+    EXPECT_REL_NEAR_F64(cephes::i1e(-0.0), 0.0);
+    // EXPECT_REL_NEAR_F64(cephes::i1e(-10.0), );
+
+    // // x <= 8.0
+    // EXPECT_REL_NEAR_F64(cephes::i1e(1.0), );
+    // EXPECT_REL_NEAR_F64(cephes::i1e(8.0), );
+
+    // // x > 8.0
+    // EXPECT_REL_NEAR_F64(cephes::i1e(9.0),  );
+    // EXPECT_REL_NEAR_F64(cephes::i1e(10.0), );
+    // EXPECT_REL_NEAR_F64(cephes::i1e(100.0), );
+}
diff --git a/tests/bessel/iv.cpp b/tests/bessel/iv.cpp
new file mode 100644
index 0000000..a83f8e3
--- /dev/null
+++ b/tests/bessel/iv.cpp
@@ -0,0 +1,26 @@
+#include <xtest.hpp>
+#include <cephes/bessel.h>
+
+/**
+    Table[NumberForm[Exp[-Abs[x]]*BesselI[nv, x], 16],
+        {nv, {0, 1}},
+        {x, {0.0, 1.0, 8.0, 9.0, 1.0, 100.}}]
+ */
+TEST(BesselIv, Branches) {
+    EXPECT_REL_NEAR_F64(cephes::iv(0, 0.0), 1.0);
+    EXPECT_REL_NEAR_F64(cephes::iv(1, 0.0), 0.0);
+    EXPECT_REL_NEAR_F64(cephes::iv(-1, 0.0), 0.0);
+
+    // x < 0
+    EXPECT_REL_NEAR_F64(cephes::iv(0, -0.0), 1.0);
+    // EXPECT_REL_NEAR_F64(cephes::iv(-10.0), );
+
+    // // x <= 8.0
+    // EXPECT_REL_NEAR_F64(cephes::iv(1.0), );
+    // EXPECT_REL_NEAR_F64(cephes::iv(8.0), );
+
+    // // x > 8.0
+    // EXPECT_REL_NEAR_F64(cephes::iv(9.0),  );
+    // EXPECT_REL_NEAR_F64(cephes::iv(10.0), );
+    // EXPECT_REL_NEAR_F64(cephes::iv(100.0), );
+}
diff --git a/tests/bessel/j0.cpp b/tests/bessel/j0.cpp
index e656c43..58975b7 100644
--- a/tests/bessel/j0.cpp
+++ b/tests/bessel/j0.cpp
@@ -1,7 +1,6 @@
-#include <limits>
-#include <complex>
-#include <cephes/bessel.h>
 #include <xtest.hpp>
+#include <cephes/bessel.h>
+
 
 TEST(BesselJ0, BasicAssertions) {
     double x, y, y_ref;
@@ -10,41 +9,41 @@ TEST(BesselJ0, Branches) {
     double x, y, y_ref;
 
     // x <= 5.0 && x < 1.0e-5
-    EXPECT_EQ(j0(0.0), 1.0);
+    EXPECT_EQ(cephes::j0(0.0), 1.0);
 
     // 1.0e-5 <= x <= 5.0
     x = 1.0;
-    y = j0(x);
+    y = cephes::j0(x);
     y_ref = 0.7651976865579666;
     XTEST_ISAPPROX_F64(y);
     x = 2.0;
-    y = j0(x);
+    y = cephes::j0(x);
     y_ref = 0.2238907791412356;
     XTEST_ISAPPROX_F64(y);
     x = 3.0;
-    y = j0(x);
+    y = cephes::j0(x);
     y_ref = -0.2600519549019333;
     XTEST_ISAPPROX_F64(y);
     x = 4.0;
-    y = j0(x);
+    y = cephes::j0(x);
     y_ref = -0.3971498098638473;
     XTEST_ISAPPROX_F64(y);
     x = 5.0;
-    y = j0(x);
+    y = cephes::j0(x);
     y_ref = -0.1775967713143385;
     XTEST_ISAPPROX_F64(y);
 
     // x > 5.0
     x = 10.0;
-    y = j0(x);
+    y = cephes::j0(x);
     y_ref = -0.2459357644513642;
     XTEST_ISAPPROX_F64(y);
     x = 100.0;
-    y = j0(x);
+    y = cephes::j0(x);
     y_ref = 0.01998585030422333;
     XTEST_ISAPPROX_F64(y);
     x = 1000.0;
-    y = j0(x);
+    y = cephes::j0(x);
     y_ref = 0.02478668615242003;
     XTEST_ISAPPROX_F64(y);
 }
diff --git a/tests/bessel/j1.cpp b/tests/bessel/j1.cpp
index bf13534..50081fd 100644
--- a/tests/bessel/j1.cpp
+++ b/tests/bessel/j1.cpp
@@ -1,40 +1,39 @@
-#include <limits>
-#include <complex>
-#include <cephes/bessel.h>
 #include <xtest.hpp>
+#include <cephes/bessel.h>
+
 
 TEST(BesselJ1, BasicAssertions) {
-    EXPECT_EQ(j1(0.0), 0.0);
-    EXPECT_EQ(j1(-0.0), 0.0);
+    EXPECT_EQ(cephes::j1(0.0), 0.0);
+    EXPECT_EQ(cephes::j1(-0.0), 0.0);
 }
 TEST(BesselJ1, Branches) {
     double x, y, y_ref;
 
     // x <= 5.0
     x = 1.0;
-    y = j1(x);
+    y = cephes::j1(x);
     y_ref = 0.4400505857449335;
     XTEST_ISAPPROX_F64(y);
     x = 2.0;
-    y = j1(x);
+    y = cephes::j1(x);
     y_ref = 0.5767248077568736;
     XTEST_ISAPPROX_F64(y);
     x = 5.0;
-    y = j1(x);
+    y = cephes::j1(x);
     y_ref = -0.3275791375914652;
     XTEST_ISAPPROX_F64(y);
 
     // x > 5.0
     x = 10.0;
-    y = j1(x);
+    y = cephes::j1(x);
     y_ref = 0.04347274616880752;
     XTEST_ISAPPROX_F64(y);
     x = 100.0;
-    y = j1(x);
+    y = cephes::j1(x);
     y_ref = -0.07714535201411228;
     XTEST_ISAPPROX_F64(y);
     x = 1000.0;
-    y = j1(x);
+    y = cephes::j1(x);
     y_ref = 0.00472831190708902;
     XTEST_ISAPPROX_F64(y);
 }
diff --git a/tests/bessel/jn.cpp b/tests/bessel/jn.cpp
new file mode 100644
index 0000000..5c07f21
--- /dev/null
+++ b/tests/bessel/jn.cpp
@@ -0,0 +1,27 @@
+#include <xtest.hpp>
+#include <cephes/bessel.h>
+
+
+TEST(BesselJn, BasicAssertions) {
+    EXPECT_REL_NEAR_F64(cephes::jn(5, 1.0), 0.0002497577302112345);
+    EXPECT_REL_NEAR_F64(cephes::jn(10, 1.0), 2.630615123687453e-10);
+    EXPECT_REL_NEAR_F64(cephes::jn(50, 1.0), 2.906004948173273e-80);
+    EXPECT_REL_NEAR_F64(cephes::jn(100, 1.0), 8.43182878962688e-189);
+}
+TEST(BesselJ1, Branches) {
+    // n < 0
+    EXPECT_REL_NEAR_F64(cephes::jn(-1, 0.0), 0.0);
+    EXPECT_TRUE(std::isnan(cephes::jn(-2, 0.0)));
+
+    // x < 0.0
+    EXPECT_REL_NEAR_F64(cephes::jn(0, -1.0), 0.7651976865579666);
+    EXPECT_REL_NEAR_F64(cephes::jn(1, -1.0), -0.4400505857449335);
+
+    // n = 0,1,2
+    EXPECT_REL_NEAR_F64(cephes::jn(0, 1.0), 0.7651976865579666);
+    EXPECT_REL_NEAR_F64(cephes::jn(1, 1.0), 0.4400505857449335);
+    EXPECT_REL_NEAR_F64(cephes::jn(2, 1.0), 0.1149034849319005);
+
+    // x < MACHEP
+    EXPECT_REL_NEAR_F64(cephes::jn(5, 0.0), 0.0);
+}
diff --git a/tests/bessel/jv.cpp b/tests/bessel/jv.cpp
new file mode 100644
index 0000000..a77dd69
--- /dev/null
+++ b/tests/bessel/jv.cpp
@@ -0,0 +1,57 @@
+#include <xtest.hpp>
+#include <cephes/bessel.h>
+
+
+/** Wolframe
+    Table[NumberForm[BesselJ[nv + 1/4., x], 16],
+        {nv, {0, 5, 10, 50, 100}},
+        {x, {1, 50}}]
+*/
+TEST(BesselJv, CoSF_Table_5p13) {
+    double nv, x;
+
+    nv = 1/4;
+    x = 1.0;
+    // TODO: large error
+    // EXPECT_REL_NEAR_F64(cephes::jv(0+nv, x),    0.7522313333407901);
+    // EXPECT_REL_NEAR_F64(cephes::jv(5+nv, x),    0.0001365611922260019);
+    // EXPECT_REL_NEAR_F64(cephes::jv(10+nv, x),   1.225744639746669e-10);
+    // EXPECT_REL_NEAR_F64(cephes::jv(50+nv, x),   9.1612835058319e-81);
+    // EXPECT_REL_NEAR_F64(cephes::jv(100+nv, x),  2.238669811158841e-189);
+
+    // nv = 1/4;
+    // x = 50.0;
+    // EXPECT_REL_NEAR_F64(cephes::jv(0+nv, x),    0.01410606268088954);
+    // EXPECT_REL_NEAR_F64(cephes::jv(5+nv, x),   -0.1041926237944591);
+    // EXPECT_REL_NEAR_F64(cephes::jv(10+nv, x),  -0.1054019881703558);
+    // EXPECT_REL_NEAR_F64(cephes::jv(50+nv, x),   0.1137789059221363);
+    // EXPECT_REL_NEAR_F64(cephes::jv(100+nv, x),  8.01262704186791e-22);
+}
+/** Wolframe
+    Table[NumberForm[BesselJ[nv + 3/4., x], 16],
+        {nv, {0, 5, 10, 50, 100}},
+        {x, {1, 50}}]
+*/
+TEST(BesselJv, CoSF_Table_5p17) {
+    double nv, x;
+
+    nv = 3/4;
+    x = 1.0;
+    // TODO: large error
+    // EXPECT_REL_NEAR_F64(cephes::jv(0+nv, x),    0.5586524932048919);
+    // EXPECT_REL_NEAR_F64(cephes::jv(5+nv, x),    0.00003952258729143145);
+    // EXPECT_REL_NEAR_F64(cephes::jv(10+nv, x),   2.615437708142973e-11);
+    // EXPECT_REL_NEAR_F64(cephes::jv(50+nv, x),   9.07139545679838e-82);
+    // EXPECT_REL_NEAR_F64(cephes::jv(100+nv, x),  1.575139282276652e-190);
+
+    // nv = 3/4;
+    // x = 50.0;
+    // EXPECT_REL_NEAR_F64(cephes::jv(0+nv, x),   -0.06874351931088636);
+    // EXPECT_REL_NEAR_F64(cephes::jv(5+nv, x),   -0.1071110937351616);
+    // EXPECT_REL_NEAR_F64(cephes::jv(10+nv, x),  -0.05465844413942518);
+    // EXPECT_REL_NEAR_F64(cephes::jv(50+nv, x),   0.0988291994940441);
+    // EXPECT_REL_NEAR_F64(cephes::jv(100+nv, x),  4.122166908740485e-22);
+}
+TEST(BesselJv, Branches) {
+
+}
diff --git a/tests/bessel/k0.cpp b/tests/bessel/k0.cpp
new file mode 100644
index 0000000..aed4e0c
--- /dev/null
+++ b/tests/bessel/k0.cpp
@@ -0,0 +1,22 @@
+#include <xtest.hpp>
+#include <cephes/bessel.h>
+
+/*
+    Table[NumberForm[BesselK[0, x], 16],
+        {x, {0.0, -10., 1.0, 8.0, 9.0, 10.0, 100.}}]
+*/
+TEST(BesselK0, Branches) {
+    EXPECT_GT(cephes::k0(0.0), 1.0e308);  // +Inf
+
+    // x < 0
+    // EXPECT_REL_NEAR_F64(cephes::k0(-10.0), );
+
+    // // x <= 8.0
+    // EXPECT_REL_NEAR_F64(cephes::k0(1.0), );
+    // EXPECT_REL_NEAR_F64(cephes::k0(8.0), );
+
+    // // x > 8.0
+    // EXPECT_REL_NEAR_F64(cephes::k0(9.0),  );
+    // EXPECT_REL_NEAR_F64(cephes::k0(10.0), );
+    // EXPECT_REL_NEAR_F64(cephes::k0(100.0), );
+}
diff --git a/tests/bessel/k0e.cpp b/tests/bessel/k0e.cpp
new file mode 100644
index 0000000..dd8e837
--- /dev/null
+++ b/tests/bessel/k0e.cpp
@@ -0,0 +1,22 @@
+#include <xtest.hpp>
+#include <cephes/bessel.h>
+
+/*
+    Table[NumberForm[Exp[-Abs[x]]*BesselK[0, x], 16],
+        {x, {0.0, -10., 1.0, 8.0, 9.0, 10.0, 100.}}]
+*/
+TEST(BesselK0Exp, Branches) {
+    EXPECT_GT(cephes::k0e(0.0), 1.0e308);  // +Inf
+
+    // x < 0
+    // EXPECT_REL_NEAR_F64(cephes::k0e(-10.0), );
+
+    // // x <= 8.0
+    // EXPECT_REL_NEAR_F64(cephes::k0e(1.0), );
+    // EXPECT_REL_NEAR_F64(cephes::k0e(8.0), );
+
+    // // x > 8.0
+    // EXPECT_REL_NEAR_F64(cephes::k0e(9.0),  );
+    // EXPECT_REL_NEAR_F64(cephes::k0e(10.0), );
+    // EXPECT_REL_NEAR_F64(cephes::k0e(100.0), );
+}
diff --git a/tests/bessel/k1.cpp b/tests/bessel/k1.cpp
new file mode 100644
index 0000000..1b88f7f
--- /dev/null
+++ b/tests/bessel/k1.cpp
@@ -0,0 +1,22 @@
+#include <xtest.hpp>
+#include <cephes/bessel.h>
+
+/*
+    Table[NumberForm[BesselK[1, x], 16],
+        {x, {0.0, -10., 1.0, 8.0, 9.0, 10.0, 100.}}]
+*/
+TEST(BesselK1, Branches) {
+    EXPECT_GT(cephes::k1(0.0), 1.0e308);  // +Inf
+
+    // x < 0
+    // EXPECT_REL_NEAR_F64(cephes::k1(-10.0), );
+
+    // // x <= 8.0
+    // EXPECT_REL_NEAR_F64(cephes::k1(1.0), );
+    // EXPECT_REL_NEAR_F64(cephes::k1(8.0), );
+
+    // // x > 8.0
+    // EXPECT_REL_NEAR_F64(cephes::k1(9.0),  );
+    // EXPECT_REL_NEAR_F64(cephes::k1(10.0), );
+    // EXPECT_REL_NEAR_F64(cephes::k1(100.0), );
+}
diff --git a/tests/bessel/k1e.cpp b/tests/bessel/k1e.cpp
new file mode 100644
index 0000000..b325f63
--- /dev/null
+++ b/tests/bessel/k1e.cpp
@@ -0,0 +1,22 @@
+#include <xtest.hpp>
+#include <cephes/bessel.h>
+
+/*
+    Table[NumberForm[Exp[-Abs[x]]*BesselK[1, x], 16],
+        {x, {0.0, -10., 1.0, 8.0, 9.0, 10.0, 100.}}]
+*/
+TEST(BesselK1Exp, Branches) {
+    EXPECT_GT(cephes::k1e(0.0), 1.0e308);  // +Inf
+
+    // x < 0
+    // EXPECT_REL_NEAR_F64(cephes::k1e(-10.0), );
+
+    // // x <= 8.0
+    // EXPECT_REL_NEAR_F64(cephes::k1e(1.0), );
+    // EXPECT_REL_NEAR_F64(cephes::k1e(8.0), );
+
+    // // x > 8.0
+    // EXPECT_REL_NEAR_F64(cephes::k1e(9.0),  );
+    // EXPECT_REL_NEAR_F64(cephes::k1e(10.0), );
+    // EXPECT_REL_NEAR_F64(cephes::k1e(100.0), );
+}
diff --git a/tests/bessel/kn.cpp b/tests/bessel/kn.cpp
new file mode 100644
index 0000000..6854897
--- /dev/null
+++ b/tests/bessel/kn.cpp
@@ -0,0 +1,25 @@
+#include <xtest.hpp>
+#include <cephes/bessel.h>
+
+/*
+    Table[NumberForm[BesselK[nv, x], 16],
+        {nv, {0, 1, 10}},
+        {x, {0.0, 1.0, 8.0, 9.0, 1.0, 100.}}]
+*/
+TEST(BesselYn, Branches) {
+    EXPECT_GT(cephes::kn(0, 0.0), 1.0e308);  // +Inf
+    EXPECT_GT(cephes::kn(1, 0.0), 1.0e308);  // +Inf
+    EXPECT_GT(cephes::kn(10, 0.0), 1.0e308);  // +Inf
+
+    // x < 0
+    // EXPECT_REL_NEAR_F64(cephes::kn(-10.0), );
+
+    // // x <= 8.0
+    // EXPECT_REL_NEAR_F64(cephes::kn(1.0), );
+    // EXPECT_REL_NEAR_F64(cephes::kn(8.0), );
+
+    // // x > 8.0
+    // EXPECT_REL_NEAR_F64(cephes::kn(9.0),  );
+    // EXPECT_REL_NEAR_F64(cephes::kn(10.0), );
+    // EXPECT_REL_NEAR_F64(cephes::kn(100.0), );
+}
diff --git a/tests/bessel/psi.cpp b/tests/bessel/psi.cpp
new file mode 100644
index 0000000..4ef34b0
--- /dev/null
+++ b/tests/bessel/psi.cpp
@@ -0,0 +1,11 @@
+#include <xtest.hpp>
+#include <cephes/bessel.h>
+
+/*
+    Table[NumberForm[PolyGamma[0, x], 16],
+        {x, {0.0, -10., 1.0, 8.0, 9.0, 10.0, 100.}}]
+*/
+TEST(DiGammaPsi, Branches) {
+    EXPECT_GT(cephes::psi(0.0), 1.0e308);  // +Inf
+
+}
diff --git a/tests/bessel/struve.cpp b/tests/bessel/struve.cpp
new file mode 100644
index 0000000..7fd215a
--- /dev/null
+++ b/tests/bessel/struve.cpp
@@ -0,0 +1,11 @@
+#include <xtest.hpp>
+#include <cephes/bessel.h>
+
+/*
+    Table[NumberForm[PolyGamma[0, x], 16],
+        {x, {0.0, -10., 1.0, 8.0, 9.0, 10.0, 100.}}]
+*/
+TEST(StruveH, Branches) {
+    // EXPECT_REL_NEAR_F64(cephes::struve(0.0, 0.0), 0.0);
+    EXPECT_REL_NEAR_F64(cephes::struve(1.0, 0.0), 0.0);
+}
diff --git a/tests/bessel/y0.cpp b/tests/bessel/y0.cpp
new file mode 100644
index 0000000..dc5f267
--- /dev/null
+++ b/tests/bessel/y0.cpp
@@ -0,0 +1,22 @@
+#include <xtest.hpp>
+#include <cephes/bessel.h>
+
+/*
+    Table[NumberForm[BesselY[0, x], 16],
+        {x, {0.0, -10., 1.0, 8.0, 9.0, 10.0, 100.}}]
+*/
+TEST(BesselY0, Branches) {
+    EXPECT_LE(cephes::y0(0.0), -1.0e308);  // -Inf
+
+    // x < 0
+    // EXPECT_REL_NEAR_F64(cephes::y0(-10.0), );
+
+    // // x <= 8.0
+    // EXPECT_REL_NEAR_F64(cephes::y0(1.0), );
+    // EXPECT_REL_NEAR_F64(cephes::y0(8.0), );
+
+    // // x > 8.0
+    // EXPECT_REL_NEAR_F64(cephes::y0(9.0),  );
+    // EXPECT_REL_NEAR_F64(cephes::y0(10.0), );
+    // EXPECT_REL_NEAR_F64(cephes::y0(100.0), );
+}
diff --git a/tests/bessel/y1.cpp b/tests/bessel/y1.cpp
new file mode 100644
index 0000000..8764f3c
--- /dev/null
+++ b/tests/bessel/y1.cpp
@@ -0,0 +1,22 @@
+#include <xtest.hpp>
+#include <cephes/bessel.h>
+
+/*
+    Table[NumberForm[BesselY[1, x], 16],
+        {x, {0.0, -10., 1.0, 8.0, 9.0, 10.0, 100.}}]
+*/
+TEST(BesselY1, Branches) {
+    EXPECT_LE(cephes::y1(0.0), -1.0e308);  // -Inf
+
+    // x < 0
+    // EXPECT_REL_NEAR_F64(cephes::y1(-10.0), );
+
+    // // x <= 8.0
+    // EXPECT_REL_NEAR_F64(cephes::y1(1.0), );
+    // EXPECT_REL_NEAR_F64(cephes::y1(8.0), );
+
+    // // x > 8.0
+    // EXPECT_REL_NEAR_F64(cephes::y1(9.0),  );
+    // EXPECT_REL_NEAR_F64(cephes::y1(10.0), );
+    // EXPECT_REL_NEAR_F64(cephes::y1(100.0), );
+}
diff --git a/tests/bessel/yn.cpp b/tests/bessel/yn.cpp
new file mode 100644
index 0000000..ed27d12
--- /dev/null
+++ b/tests/bessel/yn.cpp
@@ -0,0 +1,25 @@
+#include <xtest.hpp>
+#include <cephes/bessel.h>
+
+/*
+    Table[NumberForm[Exp[-Abs[x]]*BesselY[nv, x], 16],
+        {nv, {0, 1, 10}},
+        {x, {0.0, 1.0, 8.0, 9.0, 1.0, 100.}}]
+*/
+TEST(BesselYn, Branches) {
+    EXPECT_LE(cephes::yn(0, 0.0), -1.0e308);  // -Inf
+    EXPECT_LE(cephes::yn(1, 0.0), -1.0e308);  // -Inf
+    EXPECT_LE(cephes::yn(10, 0.0), -1.0e308);  // -Inf
+
+    // x < 0
+    // EXPECT_REL_NEAR_F64(cephes::yn(-10.0), );
+
+    // // x <= 8.0
+    // EXPECT_REL_NEAR_F64(cephes::yn(1.0), );
+    // EXPECT_REL_NEAR_F64(cephes::yn(8.0), );
+
+    // // x > 8.0
+    // EXPECT_REL_NEAR_F64(cephes::yn(9.0),  );
+    // EXPECT_REL_NEAR_F64(cephes::yn(10.0), );
+    // EXPECT_REL_NEAR_F64(cephes::yn(100.0), );
+}
diff --git a/tests/exp_int.cpp b/tests/exp_int.cpp
index 279ed10..1579345 100644
--- a/tests/exp_int.cpp
+++ b/tests/exp_int.cpp
@@ -1,35 +1,34 @@
-#include <limits>
-#include <complex>
-#include <gtest/gtest.h>
+#include <xtest.hpp>
 #include <cephes/exp_int.h>
 
+
 TEST(ExpN, Errors) {
     // n < 0
-    EXPECT_TRUE(expn(-1, 10.0) > 1e308);
+    EXPECT_TRUE(cephes::expn(-1, 10.0) > 1e308);
     // x < 0
-    EXPECT_TRUE(expn(1, -1.0) > 1e308);
+    EXPECT_TRUE(cephes::expn(1, -1.0) > 1e308);
     // x > MAXLOG ≈ 709.78
-    EXPECT_EQ(expn(1, 800), 0.0);
+    EXPECT_EQ(cephes::expn(1, 800), 0.0);
     // x==0 && n < 2
-    EXPECT_TRUE(expn(1, 0.0) > 1e308);
-    EXPECT_TRUE(expn(0, 0.0) > 1e308);
+    EXPECT_TRUE(cephes::expn(1, 0.0) > 1e308);
+    EXPECT_TRUE(cephes::expn(0, 0.0) > 1e308);
 }
 TEST(ExpN, CodecovTodo) {
     const double nan64 = std::numeric_limits<double>::quiet_NaN();
 
     // x==0.0 && n >= 2
-    EXPECT_NE(expn(2, 0.0), nan64);
-    EXPECT_NE(expn(10, 0.0), nan64);
+    EXPECT_NE(cephes::expn(2, 0.0), nan64);
+    EXPECT_NE(cephes::expn(10, 0.0), nan64);
     // n==0
-    EXPECT_NE(expn(0, 10.0), nan64);
+    EXPECT_NE(cephes::expn(0, 10.0), nan64);
     // n > 5000
-    EXPECT_NE(expn(5500, 10.0), nan64);
+    EXPECT_NE(cephes::expn(5500, 10.0), nan64);
     // n <= 5000 && x > 1.0
     // cfrac: continued fraction
-    EXPECT_NE(expn(10, 10.0), nan64);
+    EXPECT_NE(cephes::expn(10, 10.0), nan64);
     // Power series expansion
-    EXPECT_NE(expn(1, 0.5), nan64);
-    EXPECT_NE(expn(10, 0.5), nan64);
+    EXPECT_NE(cephes::expn(1, 0.5), nan64);
+    EXPECT_NE(cephes::expn(10, 0.5), nan64);
 }
 
 
@@ -38,7 +37,7 @@ TEST(SiCi, Errors) {
     double x, si, ci;
 
     // x == 0.0
-    ret = sici(0.0, &si, &ci);
+    ret = cephes::sici(0.0, &si, &ci);
     EXPECT_EQ(ret, 0);
     EXPECT_EQ(si, 0.0);
     EXPECT_LT(ci, -1e308);
@@ -49,32 +48,32 @@ TEST(SiCi, CodecovTodo) {
     double x, si, ci;
 
     // x < 0.0
-    ret = sici(-1.0, &si, &ci);
+    ret = cephes::sici(-1.0, &si, &ci);
     EXPECT_EQ(ret, 0);
     EXPECT_NE(si, nan64);
     EXPECT_NE(ci, nan64);
 
     // x > 1.0e9
-    ret = sici(1e10, &si, &ci);
+    ret = cephes::sici(1e10, &si, &ci);
     EXPECT_EQ(ret, 0);
     EXPECT_NE(si, nan64);
     EXPECT_NE(ci, nan64);
 
     // x > 4.0
     // asympt
-    ret = sici(5.0, &si, &ci);
+    ret = cephes::sici(5.0, &si, &ci);
     EXPECT_EQ(ret, 0);
     EXPECT_NE(si, nan64);
     EXPECT_NE(ci, nan64);
     // x >= 8.0
-    ret = sici(8.0, &si, &ci);
+    ret = cephes::sici(8.0, &si, &ci);
     EXPECT_EQ(ret, 0);
     EXPECT_NE(si, nan64);
     EXPECT_NE(ci, nan64);
 
     // x <= 4.0
     // asympt
-    ret = sici(3.0, &si, &ci);
+    ret = cephes::sici(3.0, &si, &ci);
     EXPECT_EQ(ret, 0);
     EXPECT_NE(si, nan64);
     EXPECT_NE(ci, nan64);
@@ -86,7 +85,7 @@ TEST(ShiChi, Errors) {
     double x, si, ci;
 
     // x == 0.0
-    ret = shichi(0.0, &si, &ci);
+    ret = cephes::shichi(0.0, &si, &ci);
     EXPECT_EQ(ret, 0);
     EXPECT_EQ(si, 0.0);
     EXPECT_LT(ci, -1e308);
@@ -97,40 +96,40 @@ TEST(ShiChi, CodecovTodo) {
     double x, si, ci;
 
     // x < 0.0
-    ret = shichi(-1.0, &si, &ci);
+    ret = cephes::shichi(-1.0, &si, &ci);
     EXPECT_EQ(ret, 0);
     EXPECT_NE(si, nan64);
     EXPECT_NE(ci, nan64);
 
     // x >= 8.0
     // chb:
-    ret = shichi(8.0, &si, &ci);
+    ret = cephes::shichi(8.0, &si, &ci);
     EXPECT_EQ(ret, 0);
     EXPECT_NE(si, nan64);
     EXPECT_NE(ci, nan64);
     // x < 18.0
-    ret = shichi(15.0, &si, &ci);
+    ret = cephes::shichi(15.0, &si, &ci);
     EXPECT_EQ(ret, 0);
     EXPECT_NE(si, nan64);
     EXPECT_NE(ci, nan64);
     // 18.0 <= x <= 88.0
-    ret = shichi(18.0, &si, &ci);
+    ret = cephes::shichi(18.0, &si, &ci);
     EXPECT_EQ(ret, 0);
     EXPECT_NE(si, nan64);
     EXPECT_NE(ci, nan64);
-    ret = shichi(88.0, &si, &ci);
+    ret = cephes::shichi(88.0, &si, &ci);
     EXPECT_EQ(ret, 0);
     EXPECT_NE(si, nan64);
     EXPECT_NE(ci, nan64);
     // x > 88.0
-    ret = shichi(100.0, &si, &ci);
+    ret = cephes::shichi(100.0, &si, &ci);
     EXPECT_EQ(ret, 0);
     EXPECT_NE(si, nan64);
     EXPECT_NE(ci, nan64);
     
     // x < 8.0
     // power series expansion
-    ret = shichi(3.0, &si, &ci);
+    ret = cephes::shichi(3.0, &si, &ci);
     EXPECT_EQ(ret, 0);
     EXPECT_NE(si, nan64);
     EXPECT_NE(ci, nan64);
diff --git a/tests/gamma/CMakeLists.txt b/tests/gamma/CMakeLists.txt
new file mode 100644
index 0000000..1ed4b9a
--- /dev/null
+++ b/tests/gamma/CMakeLists.txt
@@ -0,0 +1,15 @@
+
+add_gtest(fac)
+
+add_gtest(gamma)
+add_gtest(lgam)
+add_gtest(igam)
+add_gtest(igamc)
+add_gtest(igami)
+# add_gtest(incbet)
+# add_gtest(incbi)
+
+# add_gtest(psi)
+add_gtest(rgamma)
+
+add_gtest(beta)
diff --git a/tests/gamma/beta.cpp b/tests/gamma/beta.cpp
new file mode 100644
index 0000000..e59b0df
--- /dev/null
+++ b/tests/gamma/beta.cpp
@@ -0,0 +1,9 @@
+#include <xtest.hpp>
+#include <cephes/gamma.h>
+
+
+TEST(Beta, BasicAssertions) {
+}
+TEST(Beta, Branches) {
+    EXPECT_GT(cephes::beta(0.0, 0.0), 1.0e308);  // +Inf
+}
diff --git a/tests/gamma/fac.cpp b/tests/gamma/fac.cpp
new file mode 100644
index 0000000..0bdc40f
--- /dev/null
+++ b/tests/gamma/fac.cpp
@@ -0,0 +1,11 @@
+#include <xtest.hpp>
+#include <cephes/gamma.h>
+
+
+TEST(Fac, BasicAssertions) {
+    EXPECT_REL_NEAR_F64(cephes::fac(0.0), 1.0);
+    EXPECT_REL_NEAR_F64(cephes::fac(1.0), 1.0);
+    EXPECT_REL_NEAR_F64(cephes::fac(2.0), 2.0);
+    EXPECT_REL_NEAR_F64(cephes::fac(3.0), 6.0);
+    EXPECT_REL_NEAR_F64(cephes::fac(10.0), 3628800.0);
+}
diff --git a/tests/gamma/gamma.cpp b/tests/gamma/gamma.cpp
new file mode 100644
index 0000000..ed21d31
--- /dev/null
+++ b/tests/gamma/gamma.cpp
@@ -0,0 +1,13 @@
+#include <xtest.hpp>
+#include <cephes/gamma.h>
+
+
+TEST(Gamma, BasicAssertions) {
+    EXPECT_TRUE(std::isnan(cephes::gamma(xtest::NaN64)));
+    EXPECT_TRUE(std::isnan(cephes::gamma(0.0)));
+
+    EXPECT_REL_NEAR_F64(cephes::gamma(1.0), 1.0);
+    EXPECT_REL_NEAR_F64(cephes::gamma(2.0), 1.0);
+    EXPECT_REL_NEAR_F64(cephes::gamma(3.0), 2.0);
+    EXPECT_REL_NEAR_F64(cephes::gamma(10.0), 362880.0);
+}
diff --git a/tests/gamma/igam.cpp b/tests/gamma/igam.cpp
new file mode 100644
index 0000000..a40ec0f
--- /dev/null
+++ b/tests/gamma/igam.cpp
@@ -0,0 +1,16 @@
+#include <xtest.hpp>
+#include <cephes/gamma.h>
+
+
+TEST(GammaInc, BasicAssertions) {
+
+}
+TEST(GammaInc, Branches) {
+    // x == 0
+    EXPECT_REL_NEAR_F64(cephes::igam(1.0, 0.0), 0.0);
+    // x < 0
+    EXPECT_TRUE(std::isnan(cephes::igam(1.0, -1.0)));
+    // a <= 0
+    EXPECT_TRUE(std::isnan(cephes::igam(0.0, 1.0)));
+    EXPECT_TRUE(std::isnan(cephes::igam(-1.0, 1.0)));
+}
diff --git a/tests/gamma/igamc.cpp b/tests/gamma/igamc.cpp
new file mode 100644
index 0000000..988ed40
--- /dev/null
+++ b/tests/gamma/igamc.cpp
@@ -0,0 +1,14 @@
+#include <xtest.hpp>
+#include <cephes/gamma.h>
+
+
+TEST(GammaIncc, BasicAssertions) {
+
+}
+TEST(GammaIncc, Branches) {
+    // x < 0
+    EXPECT_TRUE(std::isnan(cephes::igamc(1.0, -1.0)));
+    // a <= 0
+    EXPECT_TRUE(std::isnan(cephes::igamc(0.0, 1.0)));
+    EXPECT_TRUE(std::isnan(cephes::igamc(-1.0, 1.0)));
+}
diff --git a/tests/gamma/igami.cpp b/tests/gamma/igami.cpp
new file mode 100644
index 0000000..ad918f6
--- /dev/null
+++ b/tests/gamma/igami.cpp
@@ -0,0 +1,9 @@
+#include <xtest.hpp>
+#include <cephes/gamma.h>
+
+
+TEST(GammaInccInv, BasicAssertions) {
+    EXPECT_REL_NEAR_F64(cephes::igami(0.5, 1.0), 0.0);
+}
+TEST(GammaInccInv, Branches) {
+}
diff --git a/tests/gamma/lgam.cpp b/tests/gamma/lgam.cpp
new file mode 100644
index 0000000..a1e0d4d
--- /dev/null
+++ b/tests/gamma/lgam.cpp
@@ -0,0 +1,11 @@
+#include <xtest.hpp>
+#include <cephes/gamma.h>
+
+
+TEST(LnGamma, BasicAssertions) {
+    EXPECT_TRUE(std::isnan(cephes::lgam(xtest::NaN64)));
+    EXPECT_TRUE(std::isinf(cephes::lgam(xtest::Inf64)));
+    EXPECT_TRUE(std::isinf(cephes::lgam(0.0)));
+
+    EXPECT_REL_NEAR_F64(cephes::lgam(1.0), 0.0);
+}
diff --git a/tests/gamma/rgamma.cpp b/tests/gamma/rgamma.cpp
new file mode 100644
index 0000000..48135a4
--- /dev/null
+++ b/tests/gamma/rgamma.cpp
@@ -0,0 +1,11 @@
+#include <xtest.hpp>
+#include <cephes/gamma.h>
+
+
+TEST(RGamma, BasicAssertions) {
+    EXPECT_REL_NEAR_F64(cephes::rgamma(0.0), 0.0);
+    EXPECT_REL_NEAR_F64(cephes::rgamma(1.0), 1.0);
+    EXPECT_REL_NEAR_F64(cephes::rgamma(2.0), 1.0);
+}
+TEST(RGamma, Branches) {
+}
diff --git a/tests/include/xtest.hpp b/tests/include/xtest.hpp
index 716d2ea..1926d49 100644
--- a/tests/include/xtest.hpp
+++ b/tests/include/xtest.hpp
@@ -4,6 +4,25 @@
 #include <limits>
 #include <gtest/gtest.h>
 
+/**
+ * @param _y_expr   Your `expr` to be evalute
+ * @param _ref_expr Reference value
+ * @param _rel_tol  Rel Tolence
+ * 
+ * - When `_ref_expr == 0`, `abs_tol = _rel_tol`
+ * - else: `abs_tol = _rel_tol * abs(_ref_expr)`
+ */
+#define EXPECT_REL_NEAR_F64_(_y_expr, _ref_expr, _rel_tol) \
+    do { \
+        double __xtest_y_ref__ = (_ref_expr); \
+        double __xtest_y_level = (__xtest_y_ref__==0) ? 1.0 : std::abs(__xtest_y_ref__); \
+        double __xtest_abs_tol = (_rel_tol) * __xtest_y_level; \
+        EXPECT_NEAR((_y_expr), __xtest_y_ref__, __xtest_abs_tol) \
+            << "rel_tol = " << _rel_tol; \
+    } while (0);
+#define EXPECT_REL_NEAR_F64(_y_expr, _ref_expr) \
+    EXPECT_REL_NEAR_F64_(_y_expr, _ref_expr, xtest::RelTolF64)
+
 /**
  * Need local var: `x`, `VAR`, `VAR_ref`
  */
@@ -39,6 +58,11 @@ constexpr T rel_tol_default() {
 template <typename T, std::enable_if_t<std::is_integral_v<T>, int> = 0>
 constexpr T rel_tol_default() { return 0; };
 
+/** Default relative tolerance for `float` type */
+const double RelTolF32 = rel_tol_default<float>();
+/** Default relative tolerance for `double` type */
+const double RelTolF64 = rel_tol_default<double>();
+
 template <typename T>
 bool isapprox(T x, T y, T rel_tol=rel_tol_default<T>()) {
     assert(rel_tol >= 0);