Add `inv_sqrt2pi` and `erfcc` to GMath

ellisowen
1 parent 32da71ef
Showing 3 changed files with 29 additions and 0 deletions
include/GMath.hpp
pyext/GMath.i
src/numerics/GMath.cpp
@@ -38,6 +38,7 @@ namespace gammalib {
     const double pihalf       =  1.570796326794896619231321691639751442099;
     const double inv_pihalf   =  0.6366197723675813430755350534900574;
     const double inv_sqrt4pi  =  0.2820947917738781434740397257803862929220;
+    const double inv_sqrt2pi  =  0.3989422804014327028632180827116826549172;
     const double pi2          = pi*pi;
     const double deg2rad      =  0.0174532925199432954743717;
     const double rad2deg      = 57.295779513082322864647722;
@@ -68,6 +69,7 @@ namespace gammalib {
     // Other functions
     double gammln(const double& arg);
+    double erfcc(const double& arg);
     double modulo(const double& v1, const double& v2);
 }
@@ -37,6 +37,7 @@ namespace gammalib {
     const double pihalf       =  1.570796326794896619231321691639751442099;
     const double inv_pihalf   =  0.6366197723675813430755350534900574;
     const double inv_sqrt4pi  =  0.2820947917738781434740397257803862929220;
+    const double inv_sqrt2pi  =  0.3989422804014327028632180827116826549172;
     const double pi2          = pi*pi;
     const double deg2rad      =  0.0174532925199432954743717;
     const double rad2deg      = 57.295779513082322864647722;
@@ -67,5 +68,6 @@ namespace gammalib {
     // Other functions
     double gammln(const double& arg);
+    double erfcc(const double& arg);
     double modulo(const double& v1, const double& v2);
 }
@@ -374,6 +374,31 @@ double gammalib::gammln(const double&amp; arg) {
 /***********************************************************************//**
+ * @brief Computes complementary error function
+ *
+ * @param[in] arg Argument.
+ * @return Complementary error function.
+ *
+ * Reference: http://en.wikipedia.org/wiki/Complementary_error_function
+ ***************************************************************************/
+double gammalib::erfcc(const double& arg) {
+   /* (C) Copr. 1986-92 Numerical Recipes Software 0@.1Y.. */
+  // Copied from Gaussian.cxx in the Fermi ScienceTools software
+   double t, z, ans;
+
+   z = std::abs(arg);
+   t = 1.0/(1.0+0.5*z);
+   ans = t*std::exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
+         t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
+         t*(-0.82215223+t*0.17087277)))))))));
+   double result = arg >= 0.0 ? ans : 2.0-ans;
+
+   // Return result
+   return result;
+}
+
+
+/***********************************************************************//**
  * @brief Returns the remainder of the division \a v1/v2.
  *
  * @param[in] v1 Argument 1.