Merge pull request #2 from andrewfowlie/dev

minor tweaks and bug fixes

src/derivatives.cpp

@@ -25,17 +25,17 @@
     Found in Mathematica from 
     `-Sum[1/n^2*Limit[D[xsq*BesselK[2,Sqrt[xsq]*n],xsq],xsq->0],{n,1,Infinity}]`    
 */
-constexpr double D1_J_B_0 = pow(M_PI, 2) / 12.;
+const double D1_J_B_0 = gsl_pow_2(M_PI) / 12.;
 /*!
     First derivatives of thermal functions with respect to
     \f$y^2\f$ at \f$y^2 \to 0\f$.
     
     Found in Mathematica from 
     `-Sum[(-1)^n/n^2*Limit[D[xsq*BesselK[2, Sqrt[xsq]*n], xsq], xsq -> 0], {n, 1, Infinity}]`    
 */
-constexpr double D1_J_F_0 = -pow(M_PI, 2) / 24.;
+const double D1_J_F_0 = -gsl_pow_2(M_PI) / 24.;
 /*! Second derivative diverges for \f$y^2\to0\f$ - return infinity */
-constexpr double INF = std::numeric_limits<double>::infinity();
+const double INF = std::numeric_limits<double>::infinity();
 
 
 // Bessel function representation of derivatives with respect to \f$y^2\f$.

src/example.cpp

@@ -1,6 +1,6 @@
 /**
     @file
-    @brief Example program that links to the `thermal_funcs` library  in `C`.
+    @brief Example program that links to the `thermal_funcs` library in `C`.
     
     To build, `make example` builds an executable `./bin/example`.
     
@@ -14,7 +14,7 @@
 
 int main() {
   /**
-      This is an example of a call to the `thermal_funcs` library  in `C`.
+      This is an example of a call to the `thermal_funcs` library in `C`.
   */
   printf("J_B = %e\n", J_B_bessel(100.));
   printf("D1_J_B_bessel = %e\n", D1_J_B_bessel(100.));

src/makefile

@@ -32,7 +32,7 @@ math.exe: thermal_funcs.o zeta.o derivatives.o math.tm math.c
 
 ../bin/example: example.cpp ../lib/thermal_funcs.so
 	# $(CXX) $(CXXFLAGS) $< -lgsl -lgslcblas -I./ thermal_funcs.o zeta.o derivatives.o -o $@
-	$(CXX) $(CXXFLAGS) $< -L$(shell echo $$PWD)/../lib -l:thermal_funcs.so -lgsl -lgslcblas -I./ -o $@
+	$(CXX) $(CXXFLAGS) $< -L$(shell echo $$PWD)/../lib -Wl,-rpath=$(shell echo $$PWD)/../lib -l:thermal_funcs.so -lgsl -lgslcblas -I./ -o $@
 
 fortran_example.o: fortran_example.f
 	$(FORTRAN) $(FFLAGS) -c $<
@@ -41,7 +41,7 @@ fortran_thermal_funcs_wrap.o: fortran_thermal_funcs_wrap.c ../lib/thermal_funcs.
 	$(CC) $(CFLAGS) -L$(shell echo $$PWD)/../lib -l:thermal_funcs.so -I./ -c $<
 
 ../bin/fortran_example: fortran_thermal_funcs_wrap.o fortran_example.o
-	$(FORTRAN) $^ -L$(shell echo $$PWD)/../lib -l:thermal_funcs.so -o $@
+	$(FORTRAN) $^ -L$(shell echo $$PWD)/../lib -Wl,-rpath=$(shell echo $$PWD)/../lib -l:thermal_funcs.so -o $@
 
 fortran: ../bin/fortran_example
 example: ../bin/example

src/math.c

@@ -3,7 +3,7 @@
     @brief Mathematica interface to thermal functions.
 */
 
-//#include <wstp.h>
+#include <wstp.h>
 #include <thermal_funcs.h>
 #include <derivatives.h>
 

src/thermal_funcs.cpp

@@ -39,28 +39,28 @@ constexpr double neg_y_squared = -1.E3;
 constexpr double pos_y_squared = 1.E3;
 // 
 /*! Term in bosonic sum defined below eq. (2.13) */
-constexpr double a_b = pow(M_PI, 2) * exp(1.5 - 2. * M_EULER);
+const double a_b = gsl_pow_2(M_PI) * exp(1.5 - 2. * M_EULER);
 /*! Term in fermionic sum defined below eq. (2.13) */
-constexpr double a_f = 16. * a_b;
+const double a_f = 16. * a_b;
 /*! \f$\pi^2\f$ */
-constexpr double M_PI_POW_2 = pow(M_PI, 2);
+const double M_PI_POW_2 = gsl_pow_2(M_PI);
 /*! \f$1/\pi^2\f$ */
-constexpr double M_PI_POW_M2 = pow(M_PI, -2);
+const double M_PI_POW_M2 = gsl_sf_pow_int(M_PI, -2);
 /*! \f$\pi^4\f$ */
-constexpr double M_PI_POW_4 = pow(M_PI, 4);
+const double M_PI_POW_4 = gsl_pow_4(M_PI);
 
 /*!
     Bosonic thermal function at \f$y^2 \to 0\f$. Found in Mathematica:
     `-Sum[1/n^2 * Limit[x^2 * BesselK[2, x *n], x -> 0], {n, 1, Infinity}]`
     or from Taylor expansions.
 */
-constexpr double J_B_0 = - pow(M_PI, 4) / 45.;
+const double J_B_0 = - gsl_pow_4(M_PI) / 45.;
 /*!
     Fermionic thermal function at \f$y^2 \to 0\f$. Found in Mathematica:
     `-Sum[(-1)^n/n^2 * Limit[x^2 * BesselK[2, x *n], x -> 0], {n, 1, Infinity}]`
     or from Taylor expansions.
 */
-constexpr double J_F_0 = 7. * pow(M_PI, 4) / 360.;
+const double J_F_0 = 7. * gsl_pow_4(M_PI) / 360.;
 
 
 // Thermal functions by numerical integration
@@ -83,7 +83,7 @@ double J_integrand(double x, double y_squared, bool bosonic) {
       @param bosonic Whether bosonic (or fermionic) function required
   */
   const double sign = 1. - 2. * static_cast<double>(bosonic);
-  const double x_squared = gsl_sf_pow_int(x, 2);
+  const double x_squared = gsl_pow_2(x);
   const double r_squared = x_squared + y_squared;
   const double abs_r = sqrt(std::abs(r_squared));
   if (r_squared >= 0.) {
@@ -180,10 +180,10 @@ double *integrand_points(double y_squared, bool bosonic) {
     const int reverse_i = n_singularity - i + 1;
 
     if (bosonic) {
-      p[reverse_i] = sqrt(-gsl_sf_pow_int(2 * i, 2) * M_PI_POW_2
+      p[reverse_i] = sqrt(-gsl_pow_2(2 * i) * M_PI_POW_2
         - y_squared);
     } else {
-      p[reverse_i] = sqrt(-gsl_sf_pow_int(2 * i - 1, 2) * M_PI_POW_2
+      p[reverse_i] = sqrt(-gsl_pow_2(2 * i - 1) * M_PI_POW_2
         - y_squared);
     }
 
@@ -313,7 +313,7 @@ double J_F_quad(double y_squared, double abs_error, double rel_error,
 // Thermal functions by Taylor expansion
 
 /*! Factor that appears in terms in sum */
-const double prefactor = 2. * gsl_sf_gamma(1.5) * gsl_sf_pow_int(4., -3) /
+const double prefactor = 2. * gsl_sf_gamma(1.5) * gsl_pow_3(0.25) /
   (gsl_sf_fact(3) * M_PI_POW_2 * M_SQRTPI);
 /*! \f$\zeta(3)\f$ */
 const double zeta_3 = gsl_sf_zeta_int(3);
@@ -337,10 +337,11 @@ double gamma_sum(double y_squared, double abs_error, double rel_error,
     }
   #endif
 
-  double factor = gsl_sf_pow_int(y_squared, 3) * prefactor;
+  double factor = gsl_pow_3(y_squared) * prefactor;
   double zeta = zeta_3;
   double term = factor * zeta;
   sum += term;
+  double fraction = 1. / 512.;
 
   for (int n = 2; n <= max_n; n += 1) {
     factor *= - y_squared * n / (n + 2.) * 0.25 * M_PI_POW_M2;
@@ -349,7 +350,8 @@ double gamma_sum(double y_squared, double abs_error, double rel_error,
     if (bosonic) {
       term = zeta * factor;
     } else {
-      term = zeta * factor * (-1. + 0.125 * gsl_sf_pow_int(0.25, n + 2));
+      fraction *= 0.25;
+      term = zeta * factor * (fraction - 1.);
     }
 
     sum += term;
@@ -407,15 +409,15 @@ double J_B_taylor(double y_squared, double abs_error, double rel_error,
   double real_y_cubed;
 
   if (y_squared >= 0.) {
-    real_y_cubed = y_squared * sqrt(y_squared);
+    real_y_cubed = pow(y_squared, 1.5);
   } else {
-    real_y_cubed = std::abs(y_squared) * sqrt(std::abs(y_squared)) * M_SQRT1_2;
+    real_y_cubed = pow(std::abs(y_squared), 1.5) * M_SQRT1_2;
   }
 
   double taylor_sum = - M_PI_POW_4 / 45.
                       + M_PI_POW_2 / 12. * y_squared
                       - M_PI / 6. * real_y_cubed
-                      - gsl_sf_pow_int(y_squared, 2)
+                      - gsl_pow_2(y_squared)
                         * log(std::abs(y_squared) / a_b) / 32.;
 
   const double sum = gamma_sum(y_squared, abs_error, rel_error, max_n, true,
@@ -455,7 +457,7 @@ double J_F_taylor(double y_squared, double abs_error, double rel_error,
 
   double taylor_sum = 7. / 360. * M_PI_POW_4
                       - M_PI_POW_2 / 24. * y_squared
-                      - gsl_sf_pow_int(y_squared, 2)
+                      - gsl_pow_2(y_squared)
                         * log(std::abs(y_squared) / a_f) / 32.;
 
   const double sum = gamma_sum(y_squared, abs_error, rel_error, max_n, false,
@@ -542,7 +544,7 @@ double bessel_sum(double y_squared, double abs_error, double rel_error,
 
   for (int n = 2; n <= max_n; n += 1) {
     const double n_double = static_cast<double>(n);
-    factor *= sign * gsl_sf_pow_int((n_double - 1.) / n_double, 2);
+    factor *= sign * gsl_pow_2((n_double - 1.) / n_double);
     const double term = factor * K2(n_double * y, fast);
     sum += term;
 

src/zeta.cpp

@@ -109,7 +109,8 @@ cdouble polylog(double s, cdouble a, int N) {
     cdouble term = a;
     cdouble sum = term;
     for (int i = 2; i <= N; i += 1) {
-      term *= a * pow((i - 1) / i, 2.5);
+      const double x = (i - 1) / static_cast<double>(i);
+      term *= a * pow(x, 2.5);
       sum += term;
     }
     return sum;

src/zeta.h

@@ -11,15 +11,7 @@
 /*! A complex double */
 typedef std::complex<double> cdouble;
 
-#ifdef __cplusplus
-extern "C" {
-#endif
-
 cdouble hurwitz_zeta(double s, cdouble a, int N = 50);
 cdouble polylog(double s, cdouble a, int N = 50);
 
-#ifdef __cplusplus
-}
-#endif
-
 #endif  // _ZETA_H_

thermal_funcs/thermal_funcs.py

@@ -39,23 +39,24 @@
        3.629909, 5.324542, 6.852961, 8.274402])
 array([-4.879475, -5.299257, -4.535294, -2.666966, -0.360079, -0.472706,
        -0.498092, -0.4295  , -0.358932, -0.298332])
-array([-4.242277, -4.526832, -4.016455, -2.663948, -0.616588, -0.472706,
-       -0.498092, -0.4295  , -0.358932, -0.298332])
+array([-4.242277, -4.526832, -4.016455, -2.663948, -0.616588, -0.506298,
+       -0.513178, -0.43695 , -0.362923, -0.3006  ])
 
 >>> for method in METHODS:
-...     np.vectorize(J_B)(y_squared, method)
-array([-1.033016, -2.37644 , -3.229206, -3.458686, -2.874375, -1.659047,
-       -1.138573, -0.837874, -0.639674, -0.50042 ])
-array([-1.033015, -2.376438, -3.229206, -3.458678, -2.874366, -1.659047,
-       -1.138573, -0.837874, -0.639674, -0.50042 ])
-array([-21.416452, -15.652599, -10.752462,  -6.674717,  -3.414974,
-        -1.766015,  -2.101377,  -3.512787,  -5.883763,  -9.171826])
+...     np.vectorize(J_F)(y_squared, method)
+array([2.553786, 4.624669, 4.576   , 3.667466, 2.445411, 1.537664,
+       1.093866, 0.81725 , 0.629048, 0.494543])
+array([2.553289, 4.624714, 4.575998, 3.667452, 2.445408, 1.537664,
+       1.093866, 0.81725 , 0.629048, 0.494543])
+array([15.433261, 11.302616,  7.689656,  4.713146,  2.555124,  1.642142,
+        1.995041,  3.228617,  5.201144,  7.83237 ])
 array([8.274402, 6.852961, 5.324542, 3.629909, 1.592409, 1.592409,
        3.629909, 5.324542, 6.852961, 8.274402])
-array([-4.879475, -5.299257, -4.535294, -2.666966, -0.360079, -0.472706,
-       -0.498092, -0.4295  , -0.358932, -0.298332])
-array([-4.242277, -4.526832, -4.016455, -2.663948, -0.616588, -0.472706,
-       -0.498092, -0.4295  , -0.358932, -0.298332])
+array([4.879475, 5.299257, 4.535294, 2.666966, 0.360079, 0.472706,
+       0.498092, 0.4295  , 0.358932, 0.298332])
+array([6.417719, 6.831623, 4.92034 , 2.349209, 0.185971, 0.446744,
+       0.484677, 0.422546, 0.355115, 0.296133])
+
 
 Test their derivatives.