Catch underflows; fix some off-by-one errors. Some test values pass now.

Former-commit-id: 3add9f21ad6fee5d443fafc3a6b2c02f2f8b6524

qpms/ewald.c

@@ -81,7 +81,6 @@ qpms_ewald32_constants_t *qpms_ewald32_constants_init(const qpms_l_t lMax /*, co
       c->s1_jMaxes[y] = -1;
   }
 
-  c->s1_constfacs[0]; //WTF???
   c->s1_constfacs_base = malloc(s1_constfacs_sz * sizeof(complex double));
   size_t s1_constfacs_sz_cumsum = 0;
   for (qpms_y_t y = 0; y < c->nelem_sc; ++y) {
@@ -176,7 +175,7 @@ int ewald32_sigma_long_shiftedpoints (
   assert(commonfac > 0);
 
   // space for Gamma_pq[j]'s
-  qpms_csf_result Gamma_pq[lMax/2];
+  qpms_csf_result Gamma_pq[lMax/2+1];
 
   // CHOOSE POINT BEGIN
   for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP
@@ -190,9 +189,9 @@ int ewald32_sigma_long_shiftedpoints (
     // R-DEPENDENT BEGIN
     complex double gamma_pq = lilgamma(rbeta_pq/k);
     complex double z = csq(gamma_pq*k/(2*eta)); // Když o tom tak přemýšlím, tak tohle je vlastně vždy reálné
-    for(qpms_l_t j = 0; j < lMax/2; ++j) {
+    for(qpms_l_t j = 0; j <= lMax/2; ++j) {
       int retval = complex_gamma_inc_e(0.5-j, z, Gamma_pq+j);
-      if(retval) abort();
+      if(!(retval==0 ||retval==GSL_EUNDRFLW)) abort();
     }
     // R-DEPENDENT END
     
@@ -206,7 +205,8 @@ int ewald32_sigma_long_shiftedpoints (
         complex double e_imalpha_pq = cexp(I*m*arg_pq);
         complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c);
         double jsum_err, jsum_err_c; kahaninit(&jsum_err, &jsum_err_c); // TODO do I really need to kahan sum errors?
-        for(qpms_l_t j = 0; j < (n-abs(m))/2; ++j) {
+        assert((n-abs(m))/2 == c->s1_jMaxes[y]);
+        for(qpms_l_t j = 0; j <= c->s1_jMaxes[y]/*(n-abs(m))/2*/; ++j) { // FIXME </<= ?
           complex double summand = pow(rbeta_pq/k, n-2*j) 
             * e_imalpha_pq  * c->legendre0[gsl_sf_legendre_array_index(n,abs(m))] // This line can actually go outside j-loop
             * cpow(gamma_pq, 2*j-1) // * Gamma_pq[j] bellow (GGG) after error computation

qpms/ewald.h

@@ -1,6 +1,7 @@
 #ifndef EWALD_H
 #define EWALD_H
 #include <gsl/gsl_sf_result.h>
+#include <gsl/gsl_errno.h>
 #include <math.h> // for inlined lilgamma
 #include <complex.h>
 #include "qpms_types.h"
@@ -24,6 +25,10 @@
  */
 
 
+// Use this handler to ignore underflows of incomplete gamma.
+gsl_error_handler_t IgnoreUnderflowsGSLErrorHandler;
+
+
 /* Object holding the constant factors from [1, (4.5)] */ 
 typedef struct { 
 	qpms_l_t lMax;

qpms/ewaldgammas.c

@@ -13,6 +13,25 @@
 #define COMPLEXPART_REL_ZERO_LIMIT 1e-14
 #endif
 
+gsl_error_handler_t IgnoreUnderflowsGSLErrorHandler;
+
+void IgnoreUnderflowsGSLErrorHandler (const char * reason, 
+                    const char * file,
+                    const int line, 
+                    const int gsl_errno) {
+  if (gsl_errno == GSL_EUNDRFLW)
+    return;
+
+  gsl_stream_printf ("ERROR", file, line, reason);
+
+  fflush(stdout);
+  fprintf (stderr, "Underflow-ignoring error handler invoked.\n");
+  fflush(stderr);
+
+  abort();
+}
+
+
 // DLMF 8.7.3 (latter expression) for complex second argument
 // BTW if a is large negative, it might take a while to evaluate.
 int cx_gamma_inc_series_e(double a, complex z, qpms_csf_result * result) {
@@ -22,8 +41,10 @@ int cx_gamma_inc_series_e(double a, complex z, qpms_csf_result * result) {
     GSL_ERROR("Undefined for non-positive integer values", GSL_EDOM);
   }
   gsl_sf_result fullgamma;
-  int retval;
-  if (GSL_SUCCESS != (retval = gsl_sf_gamma_e(a, &fullgamma))){
+  int retval = gsl_sf_gamma_e(a, &fullgamma);
+  if (GSL_EUNDRFLW == retval)
+    result->err += DBL_MIN;
+  else if (GSL_SUCCESS != retval){
     result->val = NAN + NAN*I; result->err = NAN;
     return retval;
   }

tests/ewalds.c

@@ -61,6 +61,7 @@ void ewaldtest_triang_results_free(ewaldtest_triang_results *r) {
 ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p);
 
 int main() {
+  gsl_set_error_handler(IgnoreUnderflowsGSLErrorHandler);
   for (size_t i = 0; i < sizeof(paramslist)/sizeof(ewaldtest_triang_params); ++i) {
     ewaldtest_triang_params p = paramslist[i];
     ewaldtest_triang_results *r = ewaldtest_triang(p);