From 4dd3234b0a0493c2cb6012f9f488e196a7e5aec6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 9 Oct 2019 13:20:20 +0300 Subject: [PATCH] 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 | 21 ++++++++++----------- qpms/ewald.h | 40 ++++++++++++++++++++++++++++++++++------ qpms/ewaldsf.c | 2 +- 3 files changed, 45 insertions(+), 18 deletions(-) diff --git a/qpms/ewald.c b/qpms/ewald.c index 5d10f44..b606c80 100644 --- a/qpms/ewald.c +++ b/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 diff --git a/qpms/ewald.h b/qpms/ewald.h index 2a7fec6..cb2a823 100644 --- a/qpms/ewald.h +++ b/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) diff --git a/qpms/ewaldsf.c b/qpms/ewaldsf.c index b14f07d..ddf0bc1 100644 --- a/qpms/ewaldsf.c +++ b/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; }