Precompiled tabulated version of gaunt coeffs + tests.

Former-commit-id: 3fdd540ca817b1d5cc576b54baa3e8f42bd958ce

qpms/data/gaunt18baredouble.REMOVED.git-id

@@ -0,0 +1 @@
+588d34f4170323efa5ed1965363b7fa8452075f6

qpms/data/gaunt18qmaxcumsum.REMOVED.git-id

@@ -0,0 +1 @@
+91600a82980fce28d72d0b39ee3528876f582f27

qpms/data/gaunt30qmaxcumsum.REMOVED.git-id

@@ -0,0 +1 @@
+ce0cec7c6c6bf3dce392b633833e59f96c02b2a8

qpms/gaunt.h

@@ -10,4 +10,21 @@ static inline int gaunt_q_max(int m, int n, int mu, int nu) {
 
 void gaunt_xu(int m, int n, int mu, int nu, int qmax, double *v_aq, int *err);
 //int gaunt(int m, int n, int mu, int nu, double *v_aq);
+
+
+#ifdef GAUNT_PRECOMPILED
+extern const double  gaunt_table[];
+extern const int gaunt_table_lMax;
+
+// Returns given gaunt coeff. If not in range, return nan
+// TODO distinguish between invalid input and limitations given by gaunt_table_lMax
+double gaunt_table_retrieve_single(int m, int n, int mu, int nu, int q);
+
+// returns pointer to the memory where gaunt(m, n, mu, nu, q=0) is placed
+// returns NULL if invalid input or exceeding gaunt_table_lMax
+double const * gaunt_table_retrieve_allq(int m, int n, int mu, int nu);
+
+int gaunt_table_or_xu_fill(double *target, int m, int n, int mu, int nu);
+#endif //GAUNT_PRECOMPILED
+
 #endif //GAUNT_H

qpms/gaunt_precompiled.c

@@ -0,0 +1,61 @@
+#ifndef GAUNT_PRECOMPILED
+#define GAUNT_PRECOMPILED
+#endif
+#include "gaunt.h"
+#include <assert.h>
+#include <stddef.h>
+#include <math.h>
+
+const int gaunt_table_lMax = 18;
+
+const double gaunt_table[] = {
+#include "data/gaunt18baredouble"
+};
+
+const size_t gaunt_table_qmaxcumsum[] = {
+#include "data/gaunt18qmaxcumsum"
+};
+
+/* Pořadí indexů:
+ * for (n = 0; n <= lMax; n++)
+ *   for (m = -n; m <= n; m++)
+ *     for (nu = 0; nu <= lMax; nu++)
+ *       for (mu = -nu; mu <= nu; mu++)
+ *         for (q = 0; q < min(n, nu, (n + nu - |m + mu|)/2); q++)
+ */
+
+
+double const * gaunt_table_retrieve_allq(int m, int n, int mu, int nu) {
+	assert(abs(m) <= n);
+	assert(abs(mu) <= nu);
+	if (n > gaunt_table_lMax || nu > gaunt_table_lMax)
+		return NULL;
+	size_t x = n * (size_t) (n + 1) + (ptrdiff_t) m;
+	size_t xu = nu * (size_t) (nu + 1) + (ptrdiff_t) mu;
+	size_t xcount = gaunt_table_lMax * (size_t) (gaunt_table_lMax + 2) + 1;
+	size_t idx = x * xcount + xu;
+	return gaunt_table + gaunt_table_qmaxcumsum[idx];
+}
+
+
+double gaunt_table_retrieve_single(int m, int n, int mu, int nu, int q) {
+	double const * gauntptr = gaunt_table_retrieve_allq(m, n, mu, nu);
+	if (NULL == gauntptr)
+		return NAN;
+	else return gauntptr[q];
+}
+
+int gaunt_table_or_xu_fill(double *target, int m, int n, int mu, int nu) {
+	double const * gauntptr = gaunt_table_retrieve_allq(m, n, mu, nu);
+	int qmax = gaunt_q_max(m,n,mu,nu);
+	if (gauntptr) { // table hit
+		for(int q = 0; q <= qmax; ++qmax) target[q] = gauntptr[q];
+		return 0;
+	} else {
+		int err;
+		gaunt_xu(m, n, mu, nu, qmax, target, &err);
+		return err;
+	}
+}
+
+

qpms/test_gaunt.c

@@ -0,0 +1,31 @@
+#define GAUNT_PRECOMPILED
+#include "gaunt.h"
+#include <stdio.h>
+#include <math.h>
+
+const double rerrth = 1e-12;
+
+int main()
+{
+	int lMax = gaunt_table_lMax;
+	for (int n = 1; n <= lMax; n++)
+		for (int m = -n; m <= n; m++)
+			for (int nu = 1; nu <= lMax; nu++)
+				for (int mu = -nu; mu <= nu; mu++) {
+					int err;
+					int qmax = gaunt_q_max(m,n,mu,nu);
+					double gc_xu[qmax+1];
+					gaunt_xu(m,n,mu,nu,qmax,gc_xu,&err);
+					double const * gc_table = gaunt_table_retrieve_allq(m, n, mu, nu);
+					for (int q = 0; q < qmax; ++q) {
+						double rerr = (gc_xu[q] || gc_table[q]) ? 2 * (gc_xu[q] - gc_table[q]) / fabs(gc_xu[q] + gc_table[q]) : 0;
+						printf("%.5e %s %d %d %d %d %d %.16e %.16e\n",
+								rerr, (fabs(rerr) > rerrth) ? "!" : " ", m, n, mu, nu, q, gc_xu[q], gc_table[q]);
+					}
+				}
+	return 0;
+}
+
+
+
+

tests/gauntmaxcumsum.c

@@ -0,0 +1,25 @@
+#include "../qpms/gaunt.h"
+#include <stdio.h>
+#include <stdlib.h>
+
+//const int lMax = 30;
+
+int main(int argc, char **argv) 
+{
+	int lMax;
+	if (argc < 2)
+		lMax = 30;
+	else lMax = atoi(argv[1]);
+
+	printf("// assuming lMax == %d:\n", lMax);
+	size_t qmaxcumsum = 0;
+	printf("0x%zx,\n", qmaxcumsum);
+	for (int n = 0; n <= lMax; n++)
+		for (int m = -n; m <= n; m++)
+			for (int nu = 0; nu <= lMax; nu++)
+				for (int mu = -nu; mu <= nu; mu++) {
+					qmaxcumsum += gaunt_q_max(m, n, mu, nu) + 1;
+					printf("0x%zx,\n", qmaxcumsum);
+				}
+	return 0;
+}