Merge pull request #4 from Cactus-proj/dev

cephes: init bessel && gamma tests

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)

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

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

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

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)

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)

tests/airy.cpp → tests/bessel/airy.cpp

@@ -1,15 +1,14 @@
-#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();
     int ret;
     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) {

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);
+}

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);
+}

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);
+}

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), );
+}

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), );
+}

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);
 }