From c50f40a747d960004250c597818218fe73764b2f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 27 May 2020 15:52:13 +0300 Subject: [PATCH] =?UTF-8?q?Implementation=20of=20the=20=CE=94=5Fn=20factor?= =?UTF-8?q?=20as=20a=20series;=20error=20estimates.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The error estimate for the recurrence approach is buggy. Former-commit-id: f5183eaacf6a592461f07f72e04d346a42f9fca6 --- qpms/cyewaldtest.pyx | 31 ++++++++---- qpms/ewald.c | 14 ++++-- qpms/ewald.h | 22 ++++++++- qpms/ewaldsf.c | 110 ++++++++++++++++++++++++++++++++++++++----- 4 files changed, 148 insertions(+), 29 deletions(-) diff --git a/qpms/cyewaldtest.pyx b/qpms/cyewaldtest.pyx index cf93b14..5d286ff 100644 --- a/qpms/cyewaldtest.pyx +++ b/qpms/cyewaldtest.pyx @@ -3,18 +3,29 @@ from libc.stdlib cimport malloc, free, calloc import numpy as np cdef extern from "ewald.h": - void ewald3_2_sigma_long_Delta(qpms_csf_result *target, int maxn, cdouble x, cdouble z) + void ewald3_2_sigma_long_Delta(cdouble *target, double *err, int maxn, cdouble x, cdouble z) + void ewald3_2_sigma_long_Delta_series(cdouble *target, double *err, int maxn, cdouble x, cdouble z) + void ewald3_2_sigma_long_Delta_recurrent(cdouble *target, double *err, int maxn, cdouble x, cdouble z) int complex_gamma_inc_e(double a, cdouble x, int m, qpms_csf_result *result) -def e32_Delta(int maxn, cdouble x, cdouble z): - cdef qpms_csf_result *target = malloc((maxn+1)*sizeof(qpms_csf_result)) - cdef np.ndarray[cdouble, ndim=1] target_np = np.empty((maxn+1,), dtype=complex, order='C') - ewald3_2_sigma_long_Delta(target, maxn, x, z) - cdef int i - for i in range(maxn+1): - target_np[i] = target[i].val - free(target) - return target_np +def e32_Delta(int maxn, cdouble x, cdouble z, get_err=True, method='auto'): + cdef np.ndarray[double, ndim=1] err_np + cdef double[::1] err_view + cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty((maxn+1,), dtype=complex, order='C') + cdef cdouble[::1] target_view = target_np + if get_err: + err_np = np.empty((maxn+1,), order='C') + err_view = err_np + if method == 'recurrent': + ewald3_2_sigma_long_Delta_recurrent(&target_view[0], &err_view[0] if get_err else NULL, maxn, x, z) + elif method == 'series': + ewald3_2_sigma_long_Delta_series(&target_view[0], &err_view[0] if get_err else NULL, maxn, x, z) + else: + ewald3_2_sigma_long_Delta(&target_view[0], &err_view[0] if get_err else NULL, maxn, x, z) + if get_err: + return target_np, err_np + else: + return target_np def gamma_inc(double a, cdouble x, int m=0): cdef qpms_csf_result res diff --git a/qpms/ewald.c b/qpms/ewald.c index 7cf52ca..c4e861b 100644 --- a/qpms/ewald.c +++ b/qpms/ewald.c @@ -261,7 +261,8 @@ int ewald3_21_xy_sigma_long ( complex double x; // κ*γ/(2*η) complex double factor1d = 1; // the "additional" factor for the 1D case (then it is not 1) // space for Gamma_pq[j]'s - qpms_csf_result Gamma_pq[lMax/2+1]; + complex double Gamma_pq[lMax/2+1]; + double Gamma_pq_err[lMax/2+1]; // CHOOSE POINT BEGIN // TODO mayby PGen_next_sph is not the best coordinate system choice here @@ -299,10 +300,13 @@ int ewald3_21_xy_sigma_long ( complex double x2 = x*x; if(particle_shift.z == 0) { for(qpms_l_t j = 0; j <= lMax/2; ++j) { + qpms_csf_result Gam; int retval = complex_gamma_inc_e(0.5-j, x2, // we take the branch which is principal for the Re x < 0, Im x < 0 quadrant, cf. [Linton, p. 642 in the middle] - QPMS_LIKELY(creal(x2) < 0) && !signbit(cimag(x2)) ? -1 : 0, Gamma_pq+j); + QPMS_LIKELY(creal(x2) < 0) && !signbit(cimag(x2)) ? -1 : 0, &Gam); QPMS_ENSURE_SUCCESS_OR(retval, GSL_EUNDRFLW); + Gamma_pq[j] = Gam.val; + if(err) Gamma_pq_err[j] = Gam.err; } if (latdim & LAT1D) factor1d = M_SQRT1_2 * .5 * k * gamma_pq; @@ -310,7 +314,7 @@ int ewald3_21_xy_sigma_long ( QPMS_UNTESTED; // TODO check the branches/phases! complex double z = I * k * gamma_pq * particle_shift.z; - ewald3_2_sigma_long_Delta(Gamma_pq, lMax/2, x, z); + ewald3_2_sigma_long_Delta(Gamma_pq, err ? Gamma_pq_err : NULL, lMax/2, x, z); } else { QPMS_NOT_IMPLEMENTED("1D lattices in 3D space outside of the line not implemented"); } @@ -335,9 +339,9 @@ int ewald3_21_xy_sigma_long ( * c->s1_constfacs[y][j]; if(err) { // FIXME include also other errors than Gamma_pq's relative error - kahanadd(&jsum_err, &jsum_err_c, Gamma_pq[j].err * cabs(summand)); + kahanadd(&jsum_err, &jsum_err_c, Gamma_pq_err[j] * cabs(summand)); } - summand *= Gamma_pq[j].val; // GGG + summand *= Gamma_pq[j]; // GGG ckahanadd(&jsum, &jsum_c, summand); } jsum *= phasefac * factor1d; // PFC diff --git a/qpms/ewald.h b/qpms/ewald.h index 722e838..9485a2f 100644 --- a/qpms/ewald.h +++ b/qpms/ewald.h @@ -213,9 +213,27 @@ int ewald32_sr_integral(double r, double k, double n, double eta, double *result /// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum. /** \f[ \Delta_n = \int_n^\infty t^{-1/2 - n} \exp(-t + z^2/(4t))\ud t \f] * - * NOTE: The error part is not yet implemented, NANs are returned instead! + * \bug The current choice of method, based purely on the value of \a z, might be + * unsuitable especially for big values of \a maxn. + * */ -void ewald3_2_sigma_long_Delta(qpms_csf_result *target, int maxn, complex double x, complex double z); +void ewald3_2_sigma_long_Delta(complex double *target, double *target_err, int maxn, complex double x, complex double z); + +/// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum. +/** This function always uses Kambe's (corrected) recurrent formula. + * For production, use ewald3_2_sigma_long_Delta() instead. + */ +void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *target_err, int maxn, complex double x, complex double z); + +/// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum. +/** This function always uses Taylor expansion in \a z. + * For production, use ewald3_2_sigma_long_Delta() instead. + * + * \bug The error estimate seems to be wrong (too small) at least in some cases: try + * parameters maxn = 40, z = 0.5, x = -3. This might be related to the exponential growth + * of the error. + */ +void ewald3_2_sigma_long_Delta_series(complex double *target, double *target_err, int maxn, complex double x, complex double z); // General functions acc. to [2], sec. 4.6 – currently valid for 2D and 1D lattices in 3D space diff --git a/qpms/ewaldsf.c b/qpms/ewaldsf.c index 2edfcb2..4889cae 100644 --- a/qpms/ewaldsf.c +++ b/qpms/ewaldsf.c @@ -292,26 +292,112 @@ int hyperg_2F2_series(const double a, const double b, const double c, const doub // The Delta_n factor from [Kambe II], Appendix 3 // \f[ \Delta_n = \int_n^\infty t^{-1/2 - n} \exp(-t + z^2/(4t))\ud t \f] -void ewald3_2_sigma_long_Delta(qpms_csf_result *target, int maxn, complex double x, complex double z) { +void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *err, int maxn, complex double x, complex double z) { complex double expfac = cexp(-x + 0.25 * z*z / x); complex double sqrtx = csqrt(x); // TODO check carefully, which branch is needed // These are used in the first two recurrences complex double w_plus = Faddeeva_w(+z/(2*sqrtx) + I*sqrtx, 0); complex double w_minus = Faddeeva_w(-z/(2*sqrtx) + I*sqrtx, 0); QPMS_ASSERT(maxn >= 0); - if (maxn >= 0) { - target[0].val = 0.5 * M_SQRTPI * expfac * (w_minus + w_plus); - target[0].err = NAN; //TODO - } - if (maxn >= 1) { - target[1].val = I / z * M_SQRTPI * expfac * (w_minus - w_plus); - target[1].err = NAN; //TODO - } + if (maxn >= 0) + target[0] = 0.5 * M_SQRTPI * expfac * (w_minus + w_plus); + if (maxn >= 1) + target[1] = I / z * M_SQRTPI * expfac * (w_minus - w_plus); for(int n = 1; n < maxn; ++n) { // The rest via recurrence // TODO The cpow(x, 0.5 - n) might perhaps better be replaced with a recurrently computed variant - target[n+1].val = -(4 / (z*z)) * (-(0.5 - n) * target[n].val + target[n-1].val - cpow(x, 0.5 - n) * expfac); - target[n+1].err = NAN; //TODO + target[n+1] = -(4 / (z*z)) * (-(0.5 - n) * target[n] + target[n-1] - cpow(x, 0.5 - n) * expfac); } -} + if (err) { + // The error estimates for library math functions are based on + // https://www.gnu.org/software/libc/manual/html_node/Errors-in-Math-Functions.html + // and are not guaranteed to be extremely precise. + // The error estimate for Faddeeva's functions is based on the package's authors at + // http://ab-initio.mit.edu/wiki/index.php/Faddeeva_Package + // "we find that the accuracy is typically at at least 13 significant digits in both the real and imaginary parts" + // FIXME the error estimate seems might be off by several orders of magnitude (try parameters x = -3, z = 0.5, maxn=20) + double w_plus_abs = cabs(w_plus), w_minus_abs = cabs(w_minus), expfac_abs = cabs(expfac); + double w_plus_err = w_plus_abs * 1e-13, w_minus_err = w_minus_abs * 1e-13; // LPTODO argument error contrib. + double expfac_err = expfac_abs * (4 * DBL_EPSILON); // LPTODO add argument error contrib. + double z_abs = cabs(z); + double z_err = z_abs * DBL_EPSILON; + double x_abs = cabs(x); + if (maxn >= 0) + err[0] = 0.5 * M_SQRTPI * (expfac_abs * (w_minus_err + w_plus_err) + (w_minus_abs + w_plus_abs) * expfac_err); + if (maxn >= 1) + err[1] = 2 * err[0] / z_abs + cabs(target[1]) * z_err / (z_abs*z_abs); + for(int n = 1; n < maxn; ++n) { + err[n+1] = (2 * cabs(target[n+1]) / z_abs + 4 * ((0.5+n) * err[n] + err[n-1] + + pow(x_abs, 0.5 - n) * (2*DBL_EPSILON * expfac_abs + expfac_err)) // LPTODO not ideal, pow() call is an overkill + ) * z_err / (z_abs*z_abs); + } + } +} +void ewald3_2_sigma_long_Delta_series(complex double *target, double *err, int maxn, complex double x, complex double z) { + complex double w = 0.25*z*z; + double w_abs = cabs(w); + int maxk; + if (w_abs == 0) + maxk = 0; // Delta is equal to the respective incomplete Gamma functions + else { + // Estimate a suitable maximum k, using Stirling's formula, so that w**maxk / maxk! is less than DBL_EPSILON + // This implementation is quite stupid, but it is still cheap compared to the actual computation, so LPTODO better one + maxk = 1; + double log_w_abs = log(w_abs); + while (maxk * (log_w_abs - log(maxk) + 1) >= -DBL_MANT_DIG) + ++maxk; + } + // TODO asserts on maxn, maxk + + complex double *Gammas; + double *Gammas_err = NULL, *Gammas_abs = NULL; + QPMS_CRASHING_CALLOC(Gammas, maxk+maxn+1, sizeof(*Gammas)); + if(err) { + QPMS_CRASHING_CALLOC(Gammas_err, maxk+maxn+1, sizeof(*Gammas_err)); + QPMS_CRASHING_MALLOC(Gammas_abs, (maxk+maxn+1) * sizeof(*Gammas_abs)); + } + + for(int j = 0; j <= maxn+maxk; ++j) { + qpms_csf_result g; + QPMS_ENSURE_SUCCESS(complex_gamma_inc_e(0.5-j, x, 0 /* TODO branch choice */, &g)); + Gammas[j] = g.val; + if(err) { + Gammas_abs[j] = cabs(g.val); + Gammas_err[j] = g.err; + } + } + + for(int n = 0; n <= maxn; ++n) target[n] = 0; + if(err) for(int n = 0; n <= maxn; ++n) err[n] = 0; + + complex double wpowk_over_fack = 1.; + double wpowk_over_fack_abs = 1.; + for(int k = 0; k <= maxk; ++k, wpowk_over_fack *= w/k) { // TODO? Kahan sum, Horner's method? + // Also TODO? for small n, continue for higher k if possible/needed + for(int n = 0; n <= maxn; ++n) { + target[n] += Gammas[n+k] * wpowk_over_fack; + if(err) { + // DBL_EPSILON might not be the best estimate here, but... + err[n] += wpowk_over_fack_abs * Gammas_err[n+k] + DBL_EPSILON * Gammas_abs[n+k]; + wpowk_over_fack_abs *= w_abs / (k+1); + } + } + } + + // TODO add an error estimate for the k-cutoff!!! + + free(Gammas); + free(Gammas_err); + free(Gammas_abs); +} + + +void ewald3_2_sigma_long_Delta(complex double *target, double *err, int maxn, complex double x, complex double z) { + double absz = cabs(z); + if (absz < 2.) // TODO take into account also the other parameters + ewald3_2_sigma_long_Delta_series(target, err, maxn, x, z); + else + ewald3_2_sigma_long_Delta_recurrent(target, err, maxn, x, z); +} +