fix incomplete gamma for negative x; tests

Former-commit-id: 688dc830c0b9396ceea7f503783f62d760ccf601

qpms/ewaldgammas.c

@@ -65,10 +65,11 @@ int cx_gamma_inc_series_e(double a, complex z, qpms_csf_result * result) {
 
 // incomplete gamma for complex second argument
 int complex_gamma_inc_e(double a, complex double x, qpms_csf_result *result) {
-  if (0 == fabs(cimag(x)) && // x is real; just use the real fun
+  if (creal(x) >= 0 &&
-      fabs(cimag(x)) < fabs(creal(x)) * COMPLEXPART_REL_ZERO_LIMIT) {
+      (0 == fabs(cimag(x)) || // x is real positive; just use the real fun
+      fabs(cimag(x)) < fabs(creal(x)) * COMPLEXPART_REL_ZERO_LIMIT)) {
     gsl_sf_result real_gamma_inc_result;
-    int retval = gsl_sf_gamma_inc_e(a, x, &real_gamma_inc_result);
+    int retval = gsl_sf_gamma_inc_e(a, creal(x), &real_gamma_inc_result);
     result->val = real_gamma_inc_result.val;
     result->err = real_gamma_inc_result.err;
     return retval;

tests/gsl-intgamma/genvals-incgamma.py

@@ -0,0 +1,12 @@
+import random
+import math
+
+jmu = 2
+mu = 5
+
+for i in range(10000):
+    x = abs(random.expovariate(1/mu))
+    j = math.floor(abs(random.expovariate(1/jmu)))
+    print(j, x)
+    
+            

tests/gsl-intgamma/incgamma.c

@@ -0,0 +1,28 @@
+#include <gsl/gsl_sf_gamma.h>
+#include <gsl/gsl_sf_result.h>
+#include <qpms/ewald.h>
+#include <stdio.h>
+#include <math.h>
+#include <complex.h>
+//#include <gsl/gsl_integration.h>
+#include <gsl/gsl_errno.h>
+
+
+int main(int argc, char **argv) {
+  gsl_error_handler_t * old_handler=gsl_set_error_handler_off();
+  int j;
+  double x;
+  while (scanf("%d %lf", &j, &x) == 2) {
+    printf("%d %.16g", j, x);
+    qpms_csf_result res;
+    complex double argfac = 1;
+    for (int i = 0; i < 4; ++i, argfac *= I) {
+      int retval = complex_gamma_inc_e(0.5-j, argfac * x, &res);
+      printf(" | %.16g+%.16gj %.4g %d", creal(res.val), cimag(res.val), res.err, retval);
+    }
+    putchar('\n');
+  }
+  return 0;
+}
+
+