From 7093ff3add4240720965b66a58e4ce5f3e2e38a5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Mon, 20 Aug 2018 16:29:36 +0300
Subject: [PATCH] fix incomplete gamma for negative x; tests

Former-commit-id: 688dc830c0b9396ceea7f503783f62d760ccf601
---
 qpms/ewaldgammas.c                            |  7 +++--
 tests/gsl-intgamma/genvals-incgamma.py        | 12 ++++++++
 .../{genvals.py => genvals-quadratures.py}    |  0
 tests/gsl-intgamma/incgamma.c                 | 28 +++++++++++++++++++
 4 files changed, 44 insertions(+), 3 deletions(-)
 create mode 100644 tests/gsl-intgamma/genvals-incgamma.py
 rename tests/gsl-intgamma/{genvals.py => genvals-quadratures.py} (100%)
 create mode 100644 tests/gsl-intgamma/incgamma.c

diff --git a/qpms/ewaldgammas.c b/qpms/ewaldgammas.c
index 9bc186d..930160f 100644
--- a/qpms/ewaldgammas.c
+++ b/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
-      fabs(cimag(x)) < fabs(creal(x)) * COMPLEXPART_REL_ZERO_LIMIT) {
+  if (creal(x) >= 0 &&
+      (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;
diff --git a/tests/gsl-intgamma/genvals-incgamma.py b/tests/gsl-intgamma/genvals-incgamma.py
new file mode 100644
index 0000000..33580d6
--- /dev/null
+++ b/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)
+    
+            
diff --git a/tests/gsl-intgamma/genvals.py b/tests/gsl-intgamma/genvals-quadratures.py
similarity index 100%
rename from tests/gsl-intgamma/genvals.py
rename to tests/gsl-intgamma/genvals-quadratures.py
diff --git a/tests/gsl-intgamma/incgamma.c b/tests/gsl-intgamma/incgamma.c
new file mode 100644
index 0000000..155ff81
--- /dev/null
+++ b/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;
+}
+
+