diff --git a/qpms/ewald.c b/qpms/ewald.c new file mode 100644 index 0000000..2063ef2 --- /dev/null +++ b/qpms/ewald.c @@ -0,0 +1,73 @@ +#include "ewald.h" +#include +#include "indexing.h" +#include +#include "tiny_inlines.h" + +// sloppy implementation of factorial +static inline double factorial(const int n) { + assert(n >= 0); + if (n < 0) + return 0; // should not happen in the functions below. (Therefore the assert above) + else if (n <= 20) { + double fac = 1; + for (int i = 1; i <= n; ++i) + fac *= i; + return fac; + } + else + return tgamma(n + 1); // hope it's precise and that overflow does not happen +} + + + +qpms_ewald32_constants_t *qpms_ewald32_constants_init(const qpms_l_t lMax) +{ + qpms_ewald32_constants_t *c = malloc(sizeof(qpms_ewald32_constants_t)); + //if (c == NULL) return NULL; // Do I really want to do this? + c->lMax = lMax; + c->nelem = qpms_lMax2nelem(lMax); + c->s1_jMaxes = malloc(c->nelem * sizeof(qpms_l_t)); + c->s1_constfacs = malloc(c->nelem * sizeof(complex double *)); + //if (c->s1_jMaxes == NULL) return NULL; + + //determine sizes + size_t s1_constfacs_sz = 0; + for (qpms_y_t y = 0; y < c->nelem; ++y) { + qpms_l_t n; qpms_m_t m; qpms_y2mn_p(y, &m, &n); + if ((m + n) % 2 == 0) + s1_constfacs_sz += 1 + (c->s1_jMaxes[y] = (n-abs(m))/2); + else + c->s1_jMaxes[y] = -1; + } + + c->s1_constfacs[0]; + c->s1_constfacs_base = malloc(c->nelem * sizeof(complex double)); + size_t s1_constfacs_sz_cumsum = 0; + for (qpms_y_t y = 0; y < c->nelem; ++y) { + qpms_l_t n; qpms_m_t m; qpms_y2mn_p(y, &m, &n); + if ((m + n) % 2 == 0) { + c->s1_constfacs[y] = c->s1_constfacs_base + s1_constfacs_sz_cumsum; + // and here comes the actual calculation + for (qpms_l_t j = 0; j <= c->s1_jMaxes[y]; ++j){ + c->s1_constfacs[y][j] = -0.5 * ipow(n+1) * min1pow((n+m)/2) + * sqrt((2*n + 1) * factorial(n-m) * factorial(n+m)) + * min1pow(j) * pow(0.5, n-2*j) + / (factorial(j) * factorial((n-m)/2-j) * factorial((n+m)/2-j)) + * pow(0.5, 2*j-1); + } + s1_constfacs_sz_cumsum += 1 + c->s1_jMaxes[y]; + } + else + c->s1_constfacs[y] = NULL; + } + return c; +} + +void qpms_ewald32_constants_free(qpms_ewald32_constants_t *c) { + free(c->s1_constfacs); + free(c->s1_constfacs_base); + free(c->s1_jMaxes); + free(c); +} + diff --git a/qpms/ewald.h b/qpms/ewald.h index 1f4bf7b..00ba8c5 100644 --- a/qpms/ewald.h +++ b/qpms/ewald.h @@ -1,8 +1,29 @@ #ifndef EWALD_H #define EWALD_H #include -#include +#include // for inlined lilgamma #include +#include "qpms_types.h" + +/* + * Implementation of two-dimensional lattice sum in three dimensions + * according to: + * [1] C.M. Linton, I. Thompson + * Journal of Computational Physics 228 (2009) 1815–1829 + */ + +/* Object holding the constant factors from [1, (4.5)] */ +typedef struct { + qpms_l_t lMax; + qpms_y_t nelem; + qpms_l_t *s1_jMaxes; + complex double **s1_constfacs; // indices [y][j] where j is same as in [1, (4.5)] + complex double *s1_constfacs_base; // internal pointer holding the memory for the constants +} qpms_ewald32_constants_t; + +qpms_ewald32_constants_t *qpms_ewald32_constants_init(qpms_l_t lMax); +void qpms_ewald32_constants_free(qpms_ewald32_constants_t *); + typedef struct { // as gsl_sf_result, but with complex val complex double val; @@ -10,8 +31,7 @@ typedef struct { // as gsl_sf_result, but with complex val } qpms_csf_result; -// Linton&Thompson (A.9) -// TODO put into a header file as inline +// [1, (A.9)] static inline complex double lilgamma(double t) { t = fabs(t); if (t >= 1) @@ -21,10 +41,11 @@ static inline complex double lilgamma(double t) { } +// Incomplete Gamma function as series // 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 +// 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); diff --git a/qpms/tests/s1_constfacs.c b/qpms/tests/s1_constfacs.c new file mode 100644 index 0000000..fb01bd5 --- /dev/null +++ b/qpms/tests/s1_constfacs.c @@ -0,0 +1,26 @@ +#include +#include +#include +#include +#define LMAX 10 + +int main() +{ + qpms_ewald32_constants_t *c = qpms_ewald32_constants_init(LMAX); + for (qpms_l_t l = 1; l <= LMAX; ++l) + for (qpms_m_t m = -l; m <= l; ++m) { + printf("%d %d: (", l, m); + qpms_y_t y = qpms_mn2y(m, l); + if (0 == (l+m)%2) { + for (int j = 0; j <= c->s1_jMaxes[y]; ++j) + printf("%.16g %+.16gj, ", creal(c->s1_constfacs[y][j]), + cimag(c->s1_constfacs[y][j])); + } + else assert(c->s1_constfacs[y] == NULL); + printf(")\n"); + } + return 0; +} + + + diff --git a/qpms/tests/s1_constfacs.sage b/qpms/tests/s1_constfacs.sage new file mode 100644 index 0000000..00ccbd9 --- /dev/null +++ b/qpms/tests/s1_constfacs.sage @@ -0,0 +1,16 @@ +def s1_constfacs(m, n): + if (m+n) % 2 == 1: + return [] + s1consts = list() + for j in range((n-abs(m))+1): + s1consts.append( + -I**(n+1)/2 * (-1)**((n+m)/2) * sqrt((2*n+1)*factorial(n-m)*factorial(n+m)) + * (-1)**j / 2**(n-2*j) + / (factorial(j) * factorial((n-m)/2-j) * factorial((n+m)/2-j)) / 2**(2*j-1) + ) + return s1consts + +for l in range(1, 11): + for m in range (-l, l+1): + print(l, m, [N(x, prec=66) for x in s1_constfacs(m,l)]) + diff --git a/qpms/tiny_inlines.h b/qpms/tiny_inlines.h new file mode 100644 index 0000000..70a563a --- /dev/null +++ b/qpms/tiny_inlines.h @@ -0,0 +1,23 @@ +#ifndef TINY_INLINES_H +#define TINY_INLINES_H + +static inline int min1pow(int pow) { return (pow % 2) ? -1 : 1; } + +// this has shitty precision: +// static inline complex double ipow(int x) { return cpow(I, x); } + +static inline complex double ipow(int x) { + x = ((x % 4) + 4) % 4; + switch(x) { + case 0: + return 1; + case 1: + return I; + case 2: + return -1; + case 3: + return -I; + } +} + +#endif // TINY_INLINES_H diff --git a/qpms/translations.c b/qpms/translations.c index 14307bd..148a2cb 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -6,6 +6,7 @@ #include #include #include +#include "tiny_inlines.h" #include "assert_cython_workaround.h" #include "kahansum.h" #include //abort() @@ -69,14 +70,6 @@ double qpms_legendre0(int m, int n) { return pow(2,m) * sqrtpi / tgamma(.5*n - .5*m + .5) / tgamma(.5*n-.5*m); } -static inline int min1pow(int x) { - return (x % 2) ? -1 : 1; -} - -static inline complex double ipow(int x) { - return cpow(I, x); -} - // Derivative of associated Legendre polynomial at zero argument (DLMF 14.5.2) double qpms_legendreD0(int m, int n) { return -2 * qpms_legendre0(m, n);