Remove the branch cut at negative real axis in Ewald sum.

Move the cut to the negative imaginary axis instead.

Also update comments.


Former-commit-id: a70e6a4db0ac33836d672d9ee2356ea19a899a30

qpms/ewald.c

@@ -181,8 +181,10 @@ int ewald3_sigma0(complex double *result, double *err,
     const double eta, const complex double k)
 {
   qpms_csf_result gam;
-  int retval = complex_gamma_inc_e(-0.5, -csq(k/(2*eta)), 0 /*TODO*/, &gam);
-  // FIXME DO THIS CORRECTLY gam.val = conj(gam.val); // We take the other branch, cf. [Linton, p. 642 in the middle]
+  complex double z = -csq(k/(2*eta));
+  int retval = complex_gamma_inc_e(-0.5, z, 
+     // we take the branch which is principal for the Re z < 0, Im z < 0 quadrant, cf. [Linton, p. 642 in the middle]
+     QPMS_LIKELY(creal(z) < 0) && !signbit(cimag(z)) ? -1 : 0, &gam);
   if (0 != retval)
     abort();
   *result = gam.val * c->legendre0[gsl_sf_legendre_array_index(0,0)] / 2 / M_SQRTPI;
@@ -300,11 +302,9 @@ int ewald3_21_xy_sigma_long (
       gamma_pq = clilgamma(rbeta_pq/k);
       z = csq(gamma_pq*k/(2*eta)); 
       for(qpms_l_t j = 0; j <= lMax/2; ++j) {
-        // TODO COMPLEX FIXME check the branches in the old lilgamma case
-        int retval = complex_gamma_inc_e(0.5-j, z, 0 /*TODO*/, Gamma_pq+j);
-        // we take the other branch, cf. [Linton, p. 642 in the middle]: FIXME instead use the C11 CMPLX macros and fill in -O*I part to z in the line above
-        //if(creal(z) < 0) 
-        //  Gamma_pq[j].val = conj(Gamma_pq[j].val); //FIXME as noted above
+        int retval = complex_gamma_inc_e(0.5-j, z,
+          // we take the branch which is principal for the Re z < 0, Im z < 0 quadrant, cf. [Linton, p. 642 in the middle]
+          QPMS_LIKELY(creal(z) < 0) && !signbit(cimag(z)) ? -1 : 0, Gamma_pq+j);
         if(!(retval==0 || retval==GSL_EUNDRFLW)) abort();
       }
       if (latdim & LAT1D)
@@ -443,10 +443,9 @@ int ewald3_1_z_sigma_long (
     complex double gamma_pq = clilgamma(rbeta_mu/k);  // For real beta and k this is real or pure imaginary ...
     const complex double z = csq(gamma_pq*k/(2*eta));// ... so the square (this) is in fact real.
     for(qpms_l_t j = 0; j <= lMax/2; ++j) {
-      int retval = complex_gamma_inc_e(0.5-j, z, 0 /*TODO*/, Gamma_pq+j);
-      // we take the other branch, cf. [Linton, p. 642 in the middle]: FIXME instead use the C11 CMPLX macros and fill in -O*I part to z in the line above
-      //if(creal(z) < 0) 
-      //  Gamma_pq[j].val = conj(Gamma_pq[j].val); //FIXME as noted above
+      int retval = complex_gamma_inc_e(0.5-j, z,
+          // we take the branch which is principal for the Re z < 0, Im z < 0 quadrant, cf. [Linton, p. 642 in the middle]
+          QPMS_LIKELY(creal(z) < 0) && !signbit(cimag(z)) ? -1 : 0, Gamma_pq+j);
       if(!(retval==0 || retval==GSL_EUNDRFLW)) abort();
     }
     // R-DEPENDENT END

qpms/ewald.h

@@ -133,20 +133,48 @@ static inline complex double clilgamma(complex double z) {
 }
 
 /// Incomplete Gamma function as a series.
-/** DLMF 8.7.3 (latter expression) for complex second argument */
+/** DLMF 8.7.3 (latter expression) for complex second argument.
+ *
+ * The principal value is calculated. On the negative real axis
+ * (where the function has branch cut), the sign of the imaginary
+ * part is what matters (even if it is zero). Therefore one
+ * can have
+ * `cx_gamma_inc_series_e(a, z1) != cx_gamma_inc_series_e(a, z2)`
+ * even if `z1 == z2`, because `-0 == 0` according to IEEE 754.
+ * The side of the branch cut can be determined using `signbit(creal(z))`.
+ */
 int cx_gamma_inc_series_e(double a, complex z, qpms_csf_result * result);
 
 /// Incomplete Gamma function as continued fractions.
+/** 
+ * The principal value is calculated. On the negative real axis
+ * (where the function has branch cut), the sign of the imaginary
+ * part is what matters (even if it is zero). Therefore one
+ * can have
+ * `cx_gamma_inc_CF_e(a, z1) != cx_gamma_inc_CF_e(a, z2)`
+ * even if `z1 == z2`, because `-0 == 0` according to IEEE 754.
+ * The side of the branch cut can be determined using `signbit(creal(z))`.
+ */
 int cx_gamma_inc_CF_e(double a, complex z, qpms_csf_result * result);
 
 /// Incomplete gamma for complex second argument.
-/** if x is (almost) real, it just uses gsl_sf_gamma_inc_e(). */
+/** 
+ * If x is (almost) real, it just uses gsl_sf_gamma_inc_e(). 
+ * 
+ * On the negative real axis
+ * (where the function has branch cut), the sign of the imaginary
+ * part is what matters (even if it is zero). Therefore one
+ * can have
+ * `complex_gamma_inc_e(a, z1, m) != complex_gamma_inc_e(a, z2, m)`
+ * even if `z1 == z2`, because `-0 == 0` according to IEEE 754.
+ * The side of the branch cut can be determined using `signbit(creal(z))`.
+ *
+ * Another than principal branch can be selected using non-zero \a m 
+ * argument.
+ */
 int complex_gamma_inc_e(double a, complex double x,
 	/// Branch index.
-	/** If zero, the principal value is calculated; on the negative real axis
-	 * we take the limit \f[
-	 * \Gamma(a,z)=\lim_{\epsilon\to 0+}\Gamma(\mu,z-i\epsilon).
-	 * \f]
+	/** If zero, the principal value is calculated.	 
 	 * Other branches might be chosen using non-zero \a m.
 	 * In such case, the returned value corresponds to \f[
 	 * \Gamma(a,ze^{2\pi mi})=e^{2\pi mia} \Gamma(a,z) 

qpms/ewaldsf.c

@@ -195,7 +195,7 @@ int complex_gamma_inc_e(double a, complex double x, int m, qpms_csf_result *resu
     }
     complex double f = cexp(2*m*M_PI*a*I);
     result->val *= f;
-    f = 1 - f;
+    f = -f + 1;
     result->err += cabs(f) * fullgamma.err;
     result->val += f * fullgamma.val;
   }