Basic branch selection functionality for incomplete Γ.

No modifications with effects on Ewald sums done yet.


Former-commit-id: e362341713ce306482386b9e6f2a48c336fad7a1

qpms/ewald.c

@@ -181,7 +181,7 @@ 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)), &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]
   if (0 != retval)
     abort();
@@ -301,7 +301,7 @@ int ewald3_21_xy_sigma_long (
       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, Gamma_pq+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
@@ -443,7 +443,7 @@ 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, Gamma_pq+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

qpms/ewald.h

@@ -141,7 +141,23 @@ 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(). */
-int complex_gamma_inc_e(double a, complex double x, qpms_csf_result *result);
+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]
+	 * 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) 
+	 *   + (1-e^{2\pi mia}) \Gamma(a).
+	 * \f]
+	 *
+	 * If \a a is non-positive integer, the limiting value should
+	 * be used, but this is not yet implemented!
+	 */
+	int m,
+	qpms_csf_result *result);
 
 /// Exponential integral for complex second argument.
 /** If x is (almost) positive real, it just uses gsl_sf_expint_En_e(). */

qpms/ewaldsf.c

@@ -162,26 +162,44 @@ int cx_gamma_inc_CF_e(const double a, const complex double z, qpms_csf_result *r
 }
 
 
-/// Incomplete gamma function for complex second argument.
-int complex_gamma_inc_e(double a, complex double x, qpms_csf_result *result) {
+// Incomplete gamma function for complex second argument.
+int complex_gamma_inc_e(double a, complex double x, int m, qpms_csf_result *result) {
+  int retval;
   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, creal(x), &real_gamma_inc_result);
+    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;
   } else if (creal(x) >= 0 && cabs(x) > 0.5)
-    return cx_gamma_inc_CF_e(a, x, result);
-  else if (QPMS_LIKELY(a > 0 || (a % 1.0)))
-    return cx_gamma_inc_series_e(a, x, result);
+    retval = cx_gamma_inc_CF_e(a, x, result);
+  else if (QPMS_LIKELY(a > 0 || fmod(a, 1.0)))
+    retval = cx_gamma_inc_series_e(a, x, result);
   else
   /* FIXME cx_gamma_inc_series_e() probably fails for non-positive integer a.
    * This does not matter for 2D lattices in 3D space, 
    * but it might cause problems in the other cases.
    */
     QPMS_NOT_IMPLEMENTED("Incomplete Gamma function with non-positive integer a.");
+  if (m) { // Non-principal branch.
+    /* This might be sub-optimal, as Γ(a) has probably been already evaluated
+     * somewhere in the functions called above. */
+    gsl_sf_result fullgamma;
+    int retval_fg = gsl_sf_gamma_e(a, &fullgamma);
+    if (GSL_EUNDRFLW == retval_fg)
+      fullgamma.err += DBL_MIN;
+    else if (GSL_SUCCESS != retval_fg){
+      result->val = NAN + NAN*I; result->err = NAN;
+      return GSL_ERROR_SELECT_2(retval_fg, retval);
+    }
+    complex double f = cexp(2*m*M_PI*a*I);
+    result->val *= f;
+    f = 1 - f;
+    result->err += cabs(f) * fullgamma.err;
+    result->val += f * fullgamma.val;
+  }
+  return retval;
 }
 
 // Exponential integral for complex argument; !UNTESTED! and probably not needed, as I expressed everything in terms of inc. gammas anyways.
@@ -195,7 +213,7 @@ int complex_expint_n_e(int n, complex double x, qpms_csf_result *result) {
     result->err = real_expint_result.err;
     return retval;
   } else {
-    int retval = complex_gamma_inc_e(-n+1, x, result);
+    int retval = complex_gamma_inc_e(-n+1, x, 0, result);
     complex double f = cpow(x, 2*n-2);
     result->val *= f;
     result->err *= cabs(f);

qpms/qpms_c.pyx

@@ -598,10 +598,10 @@ def linton_gamma(cdouble x):
 def linton_gamma_real(double x):
     return lilgamma(x)
 
-def gamma_inc(double a, cdouble x):
+def gamma_inc(double a, cdouble x, int m = 0):
     cdef qpms_csf_result res
     with pgsl_ignore_error(15): #15 is underflow
-        complex_gamma_inc_e(a, x, &res)
+        complex_gamma_inc_e(a, x, m, &res)
     return (res.val, res.err)
 
 def gamma_inc_series(double a, cdouble x):

qpms/qpms_cdefs.pxd

@@ -572,7 +572,7 @@ cdef extern from "ewald.h":
     cdouble clilgamma(cdouble z)
     int cx_gamma_inc_series_e(double a, cdouble x, qpms_csf_result *result)
     int cx_gamma_inc_CF_e(double a, cdouble x, qpms_csf_result *result)
-    int complex_gamma_inc_e(double a, cdouble x, qpms_csf_result *result)
+    int complex_gamma_inc_e(double a, cdouble x, int m, qpms_csf_result *result)
 
     int ewald3_sigma0(cdouble *target, double *err, const qpms_ewald3_constants_t *c, double eta, cdouble wavenumber)
     int ewald3_sigma_long(cdouble *target_sigmalr_y, double *target_sigmalr_y_err, const qpms_ewald3_constants_t *c,