From f273a46a9a943f5ac52215464144142661638581 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 20 Aug 2018 15:07:22 +0300 Subject: [PATCH] Incomplete gamma functions for complex second arguments (needed in Ewald summation) Former-commit-id: 409630d01d58f8f4e69dceb3cd59af22576acc41 --- qpms/ewald.h | 31 ++++++++ qpms/ewaldgammas.c | 78 +++++++++++++++++++ tests/{gsl-int => gsl-intgamma}/genvals.py | 0 .../quad.mathematica | 0 tests/{gsl-int => gsl-intgamma}/quadratures.c | 0 .../values.in.REMOVED.git-id | 0 .../values2.in.REMOVED.git-id | 0 7 files changed, 109 insertions(+) create mode 100644 qpms/ewald.h create mode 100644 qpms/ewaldgammas.c rename tests/{gsl-int => gsl-intgamma}/genvals.py (100%) rename tests/{gsl-int => gsl-intgamma}/quad.mathematica (100%) rename tests/{gsl-int => gsl-intgamma}/quadratures.c (100%) rename tests/{gsl-int => gsl-intgamma}/values.in.REMOVED.git-id (100%) rename tests/{gsl-int => gsl-intgamma}/values2.in.REMOVED.git-id (100%) diff --git a/qpms/ewald.h b/qpms/ewald.h new file mode 100644 index 0000000..1f4bf7b --- /dev/null +++ b/qpms/ewald.h @@ -0,0 +1,31 @@ +#ifndef EWALD_H +#define EWALD_H +#include +#include +#include + +typedef struct { // as gsl_sf_result, but with complex val + complex double val; + double err; +} qpms_csf_result; + + +// Linton&Thompson (A.9) +// TODO put into a header file as inline +static inline complex double lilgamma(double t) { + t = fabs(t); + if (t >= 1) + return sqrt(t*t - 1); + else + return -I * sqrt(1 - t*t); +} + + +// DLMF 8.7.3 (latter expression) for complex second argument +int cx_gamma_inc_series_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); + +#endif //EWALD_H diff --git a/qpms/ewaldgammas.c b/qpms/ewaldgammas.c new file mode 100644 index 0000000..9bc186d --- /dev/null +++ b/qpms/ewaldgammas.c @@ -0,0 +1,78 @@ +#include "ewald.h" +#include +#include +#include "kahansum.h" +#include +#include +//#include +#include +#include +#include + +#ifndef COMPLEXPART_REL_ZERO_LIMIT +#define COMPLEXPART_REL_ZERO_LIMIT 1e-14 +#endif + +// DLMF 8.7.3 (latter expression) for complex second argument +// BTW if a is large negative, it might take a while to evaluate. +int cx_gamma_inc_series_e(double a, complex z, qpms_csf_result * result) { + if (a <= 0 && a == (int) a) { + result->val = NAN + NAN*I; + result->err = NAN; + GSL_ERROR("Undefined for non-positive integer values", GSL_EDOM); + } + gsl_sf_result fullgamma; + int retval; + if (GSL_SUCCESS != (retval = gsl_sf_gamma_e(a, &fullgamma))){ + result->val = NAN + NAN*I; result->err = NAN; + return retval; + } + + complex double sumprefac = cpow(z, a) * cexp(-z); + double sumprefac_abs = cabs(sumprefac_abs); + complex double sum, sumc; ckahaninit(&sum, &sumc); + double err, errc; kahaninit(&err, &errc); + + bool breakswitch = false; + for (int k = 0; (!breakswitch) && (a + k + 1 <= GSL_SF_GAMMA_XMAX); ++k) { + gsl_sf_result fullgamma_ak; + if (GSL_SUCCESS != (retval = gsl_sf_gamma_e(a+k+1, &fullgamma_ak))) { + result->val = NAN + NAN*I; result->err = NAN; + return retval; + } + complex double summand = - cpow(z, k) / fullgamma_ak.val; + ckahanadd(&sum, &sumc, summand); + double summanderr = fabs(fullgamma_ak.err * cabs(summand / fullgamma_ak.val)); + // TODO add also the rounding error + kahanadd(&err, &errc, summanderr); + // TODO put some smarter cutoff break here? + if (a + k >= 18 && (cabs(summand) < err || cabs(summand) < DBL_EPSILON)) + breakswitch = true; + } + sum *= sumprefac; // Not sure if not breaking the Kahan summation here + sumc *= sumprefac; + err *= sumprefac_abs; + errc *= sumprefac_abs; + ckahanadd(&sum, &sumc, 1.); + kahanadd(&err, &errc, DBL_EPSILON); + result->err = cabs(sum) * fullgamma.err + err * fullgamma.val; + result->val = sum * fullgamma.val; // + sumc*fullgamma.val??? + if (breakswitch) + return GSL_SUCCESS; + else GSL_ERROR("Overflow; the absolute value of the z argument is probably too large.", GSL_ETOL); // maybe different error code... +} + + +// incomplete gamma for complex second argument +int complex_gamma_inc_e(double a, complex double x, qpms_csf_result *result) { + if (0 == fabs(cimag(x)) && // x is real; 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, x, &real_gamma_inc_result); + result->val = real_gamma_inc_result.val; + result->err = real_gamma_inc_result.err; + return retval; + } else + return cx_gamma_inc_series_e(a, x, result); +} + diff --git a/tests/gsl-int/genvals.py b/tests/gsl-intgamma/genvals.py similarity index 100% rename from tests/gsl-int/genvals.py rename to tests/gsl-intgamma/genvals.py diff --git a/tests/gsl-int/quad.mathematica b/tests/gsl-intgamma/quad.mathematica similarity index 100% rename from tests/gsl-int/quad.mathematica rename to tests/gsl-intgamma/quad.mathematica diff --git a/tests/gsl-int/quadratures.c b/tests/gsl-intgamma/quadratures.c similarity index 100% rename from tests/gsl-int/quadratures.c rename to tests/gsl-intgamma/quadratures.c diff --git a/tests/gsl-int/values.in.REMOVED.git-id b/tests/gsl-intgamma/values.in.REMOVED.git-id similarity index 100% rename from tests/gsl-int/values.in.REMOVED.git-id rename to tests/gsl-intgamma/values.in.REMOVED.git-id diff --git a/tests/gsl-int/values2.in.REMOVED.git-id b/tests/gsl-intgamma/values2.in.REMOVED.git-id similarity index 100% rename from tests/gsl-int/values2.in.REMOVED.git-id rename to tests/gsl-intgamma/values2.in.REMOVED.git-id