Lattice sum constants fators

Former-commit-id: 0337b5e459ad57ca81c4c903f1bc4446ec7e5566

qpms/ewald.c

@@ -0,0 +1,73 @@
+#include "ewald.h"
+#include <stdlib.h>
+#include "indexing.h"
+#include <assert.h>
+#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);
+}
+

qpms/ewald.h

@@ -1,8 +1,29 @@
 #ifndef EWALD_H
 #define EWALD_H
 #include <gsl/gsl_sf_result.h>
-#include <math.h>
+#include <math.h> // for inlined lilgamma
 #include <complex.h>
+#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);
 

qpms/tests/s1_constfacs.c

@@ -0,0 +1,26 @@
+#include <qpms/ewald.h>
+#include <qpms/indexing.h>
+#include <assert.h>
+#include <stdio.h>
+#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;
+}
+
+
+

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)])
+

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

qpms/translations.c

@@ -6,6 +6,7 @@
 #include <stdbool.h>
 #include <gsl/gsl_sf_legendre.h>
 #include <gsl/gsl_sf_bessel.h>
+#include "tiny_inlines.h"
 #include "assert_cython_workaround.h"
 #include "kahansum.h"
 #include <stdlib.h> //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);