diff --git a/qpms/data/gaunt18baredouble.REMOVED.git-id b/qpms/data/gaunt18baredouble.REMOVED.git-id new file mode 100644 index 0000000..9927db1 --- /dev/null +++ b/qpms/data/gaunt18baredouble.REMOVED.git-id @@ -0,0 +1 @@ +588d34f4170323efa5ed1965363b7fa8452075f6 \ No newline at end of file diff --git a/qpms/data/gaunt18qmaxcumsum.REMOVED.git-id b/qpms/data/gaunt18qmaxcumsum.REMOVED.git-id new file mode 100644 index 0000000..cdb6188 --- /dev/null +++ b/qpms/data/gaunt18qmaxcumsum.REMOVED.git-id @@ -0,0 +1 @@ +91600a82980fce28d72d0b39ee3528876f582f27 \ No newline at end of file diff --git a/qpms/data/gaunt30qmaxcumsum.REMOVED.git-id b/qpms/data/gaunt30qmaxcumsum.REMOVED.git-id new file mode 100644 index 0000000..73280d8 --- /dev/null +++ b/qpms/data/gaunt30qmaxcumsum.REMOVED.git-id @@ -0,0 +1 @@ +ce0cec7c6c6bf3dce392b633833e59f96c02b2a8 \ No newline at end of file diff --git a/qpms/gaunt.h b/qpms/gaunt.h index 3aa72ef..809f148 100644 --- a/qpms/gaunt.h +++ b/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 diff --git a/qpms/gaunt_precompiled.c b/qpms/gaunt_precompiled.c new file mode 100644 index 0000000..3336de1 --- /dev/null +++ b/qpms/gaunt_precompiled.c @@ -0,0 +1,61 @@ +#ifndef GAUNT_PRECOMPILED +#define GAUNT_PRECOMPILED +#endif +#include "gaunt.h" +#include +#include +#include + +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; + } +} + + diff --git a/qpms/test_gaunt.c b/qpms/test_gaunt.c new file mode 100644 index 0000000..f912b5b --- /dev/null +++ b/qpms/test_gaunt.c @@ -0,0 +1,31 @@ +#define GAUNT_PRECOMPILED +#include "gaunt.h" +#include +#include + +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; +} + + + + diff --git a/tests/gauntmaxcumsum.c b/tests/gauntmaxcumsum.c new file mode 100644 index 0000000..b8b6b3e --- /dev/null +++ b/tests/gauntmaxcumsum.c @@ -0,0 +1,25 @@ +#include "../qpms/gaunt.h" +#include +#include + +//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; +}