Constant factors in general (off-plane) Ewald 2D-in-3D sum

Former-commit-id: a8224c69682b765a36988ee62e399d97cd979f2c

notes/ewald_23_z_nonzero.lyx

@@ -233,7 +233,9 @@ For
 \begin_inset Formula 
 \begin{align*}
  & =-\frac{i^{n+1}}{k^{2}\mathscr{A}}\sqrt{\pi}2\left(\left(n-m\right)/2\right)!\left(\left(n+m\right)/2\right)!\\
- & \times\sum_{\vect K_{pq}\in\Lambda^{*}}^{'}Y_{n}^{m}\left(\frac{\pi}{2},\phi_{\vect{\beta}_{pq}}\right)\sum_{j=0}^{n-\left|m\right|}\frac{\Delta_{npq}}{j!}\left(-1\right)^{j}\left(\gamma_{pq}\right)^{2j-1}\sum_{s\overset{*}{=}j}^{\min(2j,n-\left|m\right|)}\binom{j}{2j-s}\frac{\left(-\kappa z\right)^{2j-s}\left(\beta_{pq}/k\right)^{n-s}}{\left(\frac{1}{2}\left(n-m-s\right)\right)!\left(\frac{1}{2}\left(n+m-s\right)\right)!}
+ & \times\sum_{\vect K_{pq}\in\Lambda^{*}}^{'}Y_{n}^{m}\left(\frac{\pi}{2},\phi_{\vect{\beta}_{pq}}\right)\sum_{j=0}^{n-\left|m\right|}\frac{\Delta_{npq}}{j!}\left(-1\right)^{j}\left(\gamma_{pq}\right)^{2j-1}\sum_{s\overset{*}{=}j}^{\min(2j,n-\left|m\right|)}\binom{j}{2j-s}\frac{\left(-\kappa z\right)^{2j-s}\left(\beta_{pq}/k\right)^{n-s}}{\left(\frac{1}{2}\left(n-m-s\right)\right)!\left(\frac{1}{2}\left(n+m-s\right)\right)!}\\
+ & =-\frac{i^{n+1}}{k^{2}\mathscr{A}}\sqrt{\pi}2\left(\left(n-m\right)/2\right)!\left(\left(n+m\right)/2\right)!\\
+ & \times\sum_{\vect K_{pq}\in\Lambda^{*}}^{'}Y_{n}^{m}\left(\frac{\pi}{2},\phi_{\vect{\beta}_{pq}}\right)\sum_{j=0}^{n-\left|m\right|}\Delta_{npq}\left(\gamma_{pq}\right)^{2j-1}\sum_{s\overset{*}{=}j}^{\min(2j,n-\left|m\right|)}\frac{\left(-1\right)^{j}}{\left(2j-s\right)!\left(s-j\right)!}\frac{\left(-\kappa z\right)^{2j-s}\left(\beta_{pq}/k\right)^{n-s}}{\left(\frac{1}{2}\left(n-m-s\right)\right)!\left(\frac{1}{2}\left(n+m-s\right)\right)!}
 \end{align*}
 
 \end_inset

qpms/ewald.c

@@ -155,6 +155,49 @@ qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, cons
       c->s1_constfacs[y] = NULL;
   }
 
+  // -----  New generation 2D-in-3D constants ------
+  // TODO it is not necessary to treat +|m| and -|m| cases separately
+  // N.B. currently, this is only valid for EWALD32_CONSTANTS_AGNOSTIC (NOT CHECKED!)
+  c->S1_constfacs = malloc(c->nelem_sc * sizeof(complex double *));
+  //determine sizes
+  size_t S1_constfacs_sz = 0;
+  for (qpms_y_t y = 0; y < c->nelem_sc; ++y) {
+    qpms_l_t n; qpms_m_t m; qpms_y2mn_sc_p(y, &m, &n);
+    const qpms_l_t L_M = n - abs(m);
+    for(qpms_l_t j = 0; j <= L_M; ++j) { // outer sum
+      // inner sum: j <= s <= min(2*j, n - |m|), s has the same parity as n - |m|
+      for(qpms_l_t s = j + (L_M - j) % 2; 
+          (s <= 2 * j) && (s <= L_M);
+          s += 2) {
+        ++S1_constfacs_sz;
+      }
+    }
+  }
+  c->S1_constfacs_base = malloc(S1_constfacs_sz * sizeof(complex double));
+  size_t S1_constfacs_sz_cumsum = 0; // second count
+  for (qpms_y_t y = 0; y < c->nelem_sc; ++y) {
+    qpms_l_t n; qpms_m_t m; qpms_y2mn_sc_p(y, &m, &n);
+    const complex double yfactor = -2 * ipow(n+1) * M_SQRTPI  
+              * factorial((n-m)/2) * factorial((n+m)/2);
+    c->S1_constfacs[y] = c->s1_constfacs_base + S1_constfacs_sz_cumsum;
+    size_t coeffs_per_y = 0;
+    const qpms_l_t L_M = n - abs(m);
+    for(qpms_l_t j = 0; j <= L_M; ++j) { // outer sum
+      // inner sum: j <= s <= min(2*j, n - |m|), s has the same parity as n - |m|
+      for(qpms_l_t s = j + (L_M - j) % 2; 
+          (s <= 2 * j) && (s <= L_M);
+          s += 2) {
+        c->S1_constfacs_base[S1_constfacs_sz_cumsum] = 
+             yfactor * min1pow(j) 
+             / factorial(2*j-s) / factorial(s - j) 
+             / factorial_of_half(n - m - s) / factorial_of_half(n + m - s);
+        ++S1_constfacs_sz_cumsum;
+      }
+    }
+    S1_constfacs_sz_cumsum += coeffs_per_y;
+  }
+  QPMS_ASSERT(S1_constfacs_sz_cumsum = S1_constfacs_sz);
+
   // ------ the "z-axis constants" -----
   // determine sizes
   size_t s1_constfacs_1Dz_sz;

qpms/ewald.h

@@ -90,7 +90,7 @@ typedef struct qpms_ewald3_constants_t {
 	/* These are the actual numbers now: (in the EWALD32_CONSTANTS_AGNOSTIC version)
 	 * for m + n EVEN:
 	 *
-	 * s1_constfacs[y(m,n)][x(j,s)] = 
+	 * S1_constfacs[y(m,n)][x(j,s)] = 
 	 *
 	 *   -2 * I**(n+1) * sqrt(π) * ((n-m)/2)! * ((n+m)/2)! * (-1)**j /   j    \
 	 *   ----------------------------------------------------------- |        |
@@ -98,7 +98,7 @@ typedef struct qpms_ewald3_constants_t {
 	 *
 	 * for m + n ODD:
 	 *
-	 * s1_constfacs[y(m,n)][j] = 0
+	 * S1_constfacs[y(m,n)][j] = 0
 	 */
 	complex double *S1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors.
 	/// The constant factors for the long range part of a 1D Ewald sum along the \a z axis.