Hypergeometric 2F2 for calculating lattice sum errors

Former-commit-id: eac72db084c3b3bd387115255590df4269aca76f

qpms/ewald.h

@@ -76,6 +76,11 @@ int cx_gamma_inc_series_e(double a, complex z, qpms_csf_result * result);
 // if x is (almost) real, it just uses gsl_sf_gamma_inc_e
 int complex_gamma_inc_e(double a, complex double x, qpms_csf_result *result);
 
+
+// hypergeometric 2F2, used to calculate some errors
+int hyperg_2F2_series(const double a, const double b, const double c, const double d,
+    const double x, gsl_sf_result *result);
+
 #if 0
 // The integral from (4.6); maybe should be static and not here.
 int ewald32_sr_integral(double r, double k, double n, double eta, double *result, double *err, gsl_integration_workspace *workspace);

qpms/ewaldsf.c

@@ -1,6 +1,7 @@
 #include "ewald.h"
 #include <gsl/gsl_sf_gamma.h>
 #include <gsl/gsl_sf_result.h>
+#include <gsl/gsl_machine.h> // Maybe I should rather use DBL_EPSILON instead of GSL_DBL_EPSILON.
 #include "kahansum.h"
 #include <math.h>
 #include <complex.h>
@@ -98,3 +99,74 @@ int complex_gamma_inc_e(double a, complex double x, qpms_csf_result *result) {
     return cx_gamma_inc_series_e(a, x, result);
 }
 
+
+// inspired by GSL's hyperg_2F1_series
+int hyperg_2F2_series(const double a, const double b, const double c, const double d,
+    const double x, gsl_sf_result *result
+    )
+{
+  double sum_pos = 1.0;
+  double sum_neg = 0.0;
+  double del_pos = 1.0;
+  double del_neg = 0.0;
+  double del = 1.0;
+  double del_prev;
+  double k = 0.0;
+  int i = 0;
+
+  if(fabs(c) < GSL_DBL_EPSILON || fabs(d) < GSL_DBL_EPSILON) {
+    result->val = NAN; 
+    result->err = INFINITY;
+    GSL_ERROR ("error", GSL_EDOM);
+  }
+
+  do {
+    if(++i > 30000) {
+      result->val  = sum_pos - sum_neg;
+      result->err  = del_pos + del_neg;
+      result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
+      result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k)+1.0) * fabs(result->val);
+      GSL_ERROR ("error", GSL_EMAXITER);
+    }
+    del_prev = del;
+    del *= (a+k)*(b+k) * x / ((c+k) * (d+k) * (k+1.0));  /* Gauss series */
+
+    if(del > 0.0) {
+      del_pos  =  del;
+      sum_pos +=  del;
+    }
+    else if(del == 0.0) {
+      /* Exact termination (a or b was a negative integer).
+       */
+      del_pos = 0.0;
+      del_neg = 0.0;
+      break;
+    }
+    else {
+      del_neg  = -del;
+      sum_neg -=  del;
+    }
+
+    /*
+     * This stopping criteria is taken from the thesis
+     * "Computation of Hypergeometic Functions" by J. Pearson, pg. 31
+     * (http://people.maths.ox.ac.uk/porterm/research/pearson_final.pdf)
+     * and fixes bug #45926
+     */
+    if (fabs(del_prev / (sum_pos - sum_neg)) < GSL_DBL_EPSILON &&
+        fabs(del / (sum_pos - sum_neg)) < GSL_DBL_EPSILON)
+      break;
+
+    k += 1.0;
+  } while(fabs((del_pos + del_neg)/(sum_pos-sum_neg)) > GSL_DBL_EPSILON);
+
+  result->val  = sum_pos - sum_neg;
+  result->err  = del_pos + del_neg;
+  result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
+  result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k) + 1.0) * fabs(result->val);
+
+  return GSL_SUCCESS;
+}
+
+
+

qpms/tiny_inlines.h

@@ -7,6 +7,8 @@ static inline int min1pow(int pow) { return (pow % 2) ? -1 : 1; }
 
 // This is useful for calculating spherical harmonics with negative m
 // if spharm-normalised legendre functions for positive m are available.
+// TODO: write a function that gets legendre buffer, m, n, and returns the correct spharm
+// and use it in the code (mainly translations.c, ewald.c).
 static inline int min1pow_m_neg(int m) {
 	return (m < 0) ? min1pow(m) : 1;
 }

tests/gsl-intgamma/ewald2F2.c

@@ -0,0 +1,21 @@
+// c99 -I ../.. ewald2F2.c ../../qpms/ewaldsf.c -lm -lgsl -lblas
+#include <gsl/gsl_sf_result.h>
+#include <qpms/ewald.h>
+#include <stdio.h>
+#include <math.h>
+#include <gsl/gsl_errno.h>
+
+
+int main(int argc, char **argv) {
+  gsl_error_handler_t * old_handler=gsl_set_error_handler_off();
+  double a, b, c, d, x;
+  while (scanf("%lf %lf %lf %lf %lf", &a, &b, &c, &d, &x) == 5) {
+    printf("%.16g %.16g %.16g %.16g %.16g", a, b, c, d, x);
+    gsl_sf_result res;
+    int retval = hyperg_2F2_series(a, b, c, d, x, &res);
+    printf(" | %.16g (%.3g) %d\n", res.val,  res.err, retval);
+  }
+  return 0;
+}
+
+