qpms/gaunt_precompiled.c

   1 #ifndef GAUNT_PRECOMPILED
   2 #define GAUNT_PRECOMPILED
   3 #endif
   4 #include "gaunt.h"
   5 #include <assert.h>
   6 #include <stddef.h>
   7 #include <math.h>
   8
   9 const int gaunt_table_lMax = 18;
  10
  11 const double gaunt_table[] = {
  12 #include "data/gaunt18baredouble"
  13 };
  14
  15 const size_t gaunt_table_qmaxcumsum[] = {
  16 #include "data/gaunt18qmaxcumsum"
  17 };
  18
  19 /* Pořadí indexů:
  20  * for (n = 0; n <= lMax; n++)
  21  *   for (m = -n; m <= n; m++)
  22  *     for (nu = 0; nu <= lMax; nu++)
  23  *       for (mu = -nu; mu <= nu; mu++)
  24  *         for (q = 0; q < min(n, nu, (n + nu - |m + mu|)/2); q++)
  25  */
  26
  27
  28 double const * gaunt_table_retrieve_allq(int m, int n, int mu, int nu) {
  29         assert(abs(m) <= n);
  30         assert(abs(mu) <= nu);
  31         if (n > gaunt_table_lMax || nu > gaunt_table_lMax)
  32                 return NULL;
  33         size_t x = n * (size_t) (n + 1) + (ptrdiff_t) m;
  34         size_t xu = nu * (size_t) (nu + 1) + (ptrdiff_t) mu;
  35         size_t xcount = gaunt_table_lMax * (size_t) (gaunt_table_lMax + 2) + 1;
  36         size_t idx = x * xcount + xu;
  37         return gaunt_table + gaunt_table_qmaxcumsum[idx];
  38 }
  39
  40
  41 double gaunt_table_retrieve_single(int m, int n, int mu, int nu, int q) {
  42         double const * gauntptr = gaunt_table_retrieve_allq(m, n, mu, nu);
  43         if (NULL == gauntptr)
  44                 return NAN;
  45         else return gauntptr[q];
  46 }
  47
  48 int gaunt_table_or_xu_fill(double *target, int m, int n, int mu, int nu) {
  49         double const * gauntptr = gaunt_table_retrieve_allq(m, n, mu, nu);
  50         int qmax = gaunt_q_max(m,n,mu,nu);
  51         if (gauntptr) { // table hit
  52                 for(int q = 0; q <= qmax; ++qmax) target[q] = gauntptr[q];
  53                 return 0;
  54         } else {
  55                 int err;
  56                 gaunt_xu(m, n, mu, nu, qmax, target, &err);
  57                 return err;
  58         }
  59 }
  60
  61