From e7f0e0131f191c820740d73a40c933a8ebe65eca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sun, 6 Oct 2019 01:30:29 +0300 Subject: [PATCH] =?UTF-8?q?Basic=20branch=20selection=20functionality=20fo?= =?UTF-8?q?r=20incomplete=20=CE=93.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit No modifications with effects on Ewald sums done yet. Former-commit-id: e362341713ce306482386b9e6f2a48c336fad7a1 --- qpms/ewald.c | 6 +++--- qpms/ewald.h | 18 +++++++++++++++++- qpms/ewaldsf.c | 34 ++++++++++++++++++++++++++-------- qpms/qpms_c.pyx | 4 ++-- qpms/qpms_cdefs.pxd | 2 +- 5 files changed, 49 insertions(+), 15 deletions(-) diff --git a/qpms/ewald.c b/qpms/ewald.c index 3d08de0..5d10f44 100644 --- a/qpms/ewald.c +++ b/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 diff --git a/qpms/ewald.h b/qpms/ewald.h index f6138f6..2a7fec6 100644 --- a/qpms/ewald.h +++ b/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(). */ diff --git a/qpms/ewaldsf.c b/qpms/ewaldsf.c index e0596cb..b14f07d 100644 --- a/qpms/ewaldsf.c +++ b/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); diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 70c17b6..d1acf59 100644 --- a/qpms/qpms_c.pyx +++ b/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): diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index d87d27e..8b8a4ba 100644 --- a/qpms/qpms_cdefs.pxd +++ b/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,