qpms/bessel.c

   1 #include <assert.h>
   2 #include "qpms_specfunc.h"
   3 #include <stdlib.h>
   4 #include <stddef.h>
   5 #include <string.h>
   6 #include "kahansum.h"
   7 #include <gsl/gsl_sf_bessel.h>
   8 #include <complex.h>
   9 #include "qpms_error.h"
  10 #include <camos.h>
  11 #include <math.h>
  12
  13 #ifndef M_LN2
  14 #define M_LN2   0.69314718055994530942  /* log_e 2 */
  15 #endif
  16
  17 static inline complex double ipow(int x) {
  18   return cpow(I,x);
  19 }
  20
  21 #if 0
  22 // Inspired by scipy/special/_spherical_bessel.pxd
  23 static inline complex double spherical_jn(qpms_l_t l, complex double z) {
  24   if (isnan(creal(z)) || isnan(cimag(z))) return NAN+I*NAN;
  25   if (l < 0) QPMS_WTF;
  26   if (fpclassify(creal(z)) == FP_INFINITE)
  27     if(0 == cimag(z))
  28       return 0;
  29     else
  30       return INFINITY + I * INFINITY;
  31   if (z == 0)
  32     if (l == 0) return 1;
  33     else return 0;
  34   return csqrt(M_PI_2/z) * zbesj(l + .5, z);
  35 }
  36
  37 static inline complex double spherical_yn(qpms_l_t l, complex double z) {
  38   if (isnan(creal(z)) || isnan(cimag(z))) return NAN+I*NAN;
  39   if (l < 0) QPMS_WTF;
  40   if (fpclassify(creal(z)) == FP_INFINITE)
  41     if(0 == cimag(z))
  42       return 0;
  43     else
  44       return INFINITY + I * INFINITY;
  45   if (z == 0)
  46     return NAN;
  47   return csqrt(M_PI_2/z) * zbesy(l + .5, z);
  48 }
  49 #endif
  50
  51 // There is a big issue with gsl's precision of spherical bessel function; these have to be implemented differently
  52 qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double x, complex double *result_array) {
  53   int retval;
  54   double tmparr[lmax+1];
  55   if (typ != QPMS_BESSEL_REGULAR && x == 0) {
  56     for (int l = 0; l <= lmax; ++l) result_array[l] = NAN + I*NAN;
  57     QPMS_WARN("Evaluating singular Bessel functions at zero!");
  58     return QPMS_SUCCESS; // Really?
  59   }
  60   switch(typ) {
  61     case QPMS_BESSEL_REGULAR:
  62       retval = gsl_sf_bessel_jl_steed_array(lmax, x, tmparr);
  63       for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l];
  64       return retval;
  65       break;
  66     case QPMS_BESSEL_SINGULAR: //FIXME: is this precise enough? Would it be better to do it one-by-one?
  67       retval = gsl_sf_bessel_yl_array(lmax,x,tmparr);
  68       for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l];
  69       return retval;
  70       break;
  71     case QPMS_HANKEL_PLUS:
  72     case QPMS_HANKEL_MINUS:
  73       retval = gsl_sf_bessel_jl_steed_array(lmax, x, tmparr);
  74       for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l];
  75       if(retval) return retval;
  76       retval = gsl_sf_bessel_yl_array(lmax, x, tmparr);
  77       if (typ==QPMS_HANKEL_PLUS)
  78         for (int l = 0; l <= lmax; ++l) result_array[l] += I * tmparr[l];
  79       else
  80         for (int l = 0; l <= lmax; ++l) result_array[l] +=-I * tmparr[l];
  81       return retval;
  82       break;
  83     default:
  84       abort();
  85       //return GSL_EDOM;
  86   }
  87   assert(0);
  88 }
  89
  90 // TODO DOC
  91 qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, complex double x, complex double *res) {
  92   if(!cimag(x) &&
  93         creal(x) < 10000) // For large arguments (around 30000 or more), gsl_sf_bessel_jl_steed_array fails
  94     return qpms_sph_bessel_realx_fill(typ, lmax, creal(x), res);
  95   else if (isnan(creal(x)) || isnan(cimag(x)))
  96     for(qpms_l_t l = 0; l <= lmax; ++l) res[l] = NAN + I*NAN;
  97   else if (lmax < 0) QPMS_WTF;
  98   else if (fpclassify(creal(x)) == FP_INFINITE)
  99     for(qpms_l_t l = 0; l <= lmax; ++l) res[l] = INFINITY + I * INFINITY;
 100   else {
 101 try_again: ;
 102     int retry_counter = 0;
 103     const double zr = creal(x), zi = cimag(x), fnu = 0.5;
 104     const int n = lmax + 1, kode = 1 /* No exponential scaling */;
 105     double cyr[n], cyi[n];
 106     int ierr, nz;
 107     unsigned int kindchar; // Only for error output
 108     const complex double prefac = csqrt(M_PI_2/x);
 109     switch(typ) {
 110       case QPMS_BESSEL_REGULAR:
 111         kindchar = 'j';
 112         ierr = camos_zbesj(zr, zi, fnu, kode, n, cyr, cyi, &nz);
 113         break;
 114       case QPMS_BESSEL_SINGULAR:
 115         kindchar = 'y';
 116         {
 117           double cwrkr[lmax + 1], cwrki[lmax + 1];
 118           ierr = camos_zbesy(zr, zi, fnu, kode, n, cyr, cyi, &nz, cwrkr, cwrki);
 119         }
 120         break;
 121       case QPMS_HANKEL_PLUS:
 122       case QPMS_HANKEL_MINUS:
 123         kindchar = 'h';
 124         {
 125           const int m = (typ == QPMS_HANKEL_PLUS) ? 1 : 2;
 126           ierr = camos_zbesh(zr, zi, fnu, kode, m, n, cyr, cyi, &nz);
 127         }
 128         break;
 129       default:
 130         QPMS_WTF;
 131     }
 132     // TODO check for underflows? (nz != 0)
 133     if (ierr == 0 || ierr == 3) {
 134       for (qpms_l_t l = 0; l <= lmax; ++l)
 135         res[l] = prefac * (cyr[l] + I * cyi[l]);
 136       if (ierr == 3)
 137         if (creal(x) >= 10000) {
 138           QPMS_WARN_ONCE("Due to large argument, "
 139                 "Amos's zbes%c computation done but losses of significance "
 140                 "by argument reduction produce less than half of machine accuracy. "
 141                 "This warning is shown only once.",
 142                 kindchar);
 143         } else {
 144           QPMS_WARN("Amos's zbes%c computation done but losses of significance "
 145                 "by argument reduction produce less than half of machine accuracy.",
 146                 kindchar);
 147         }
 148       return QPMS_SUCCESS; //TODO maybe something else if ierr == 3
 149     }
 150     else {
 151       if (retry_counter < 5) {
 152         QPMS_WARN("Amos's zbes%c failed with ierr == %d (lMax = %d, x = %+.16g%+.16gi). Retrying.\n",
 153            kindchar, (int) ierr, lmax, creal(x), cimag(x));
 154         ++retry_counter;
 155         goto try_again;
 156       } else QPMS_PR_ERROR("Amos's zbes%c failed with ierr == %d (lMax = %d, x = %+.16g%+.16gi).",
 157            kindchar, (int) ierr, lmax, creal(x), cimag(x));
 158     }
 159   }
 160   return QPMS_SUCCESS;
 161 }
 162
 163 static inline ptrdiff_t akn_index(qpms_l_t n, qpms_l_t k) {
 164   assert(k <= n);
 165   return ((ptrdiff_t) n + 1) * n / 2 + k;
 166 }
 167 static inline ptrdiff_t bkn_index(qpms_l_t n, qpms_l_t k) {
 168   assert(k <= n+1);
 169   return ((ptrdiff_t) n + 2) * (n + 1) / 2 - 1 + k;
 170 }
 171
 172 static inline qpms_errno_t qpms_sbessel_calculator_ensure_lMax(qpms_sbessel_calculator_t *c, qpms_l_t lMax) {
 173   if (lMax <= c->lMax)
 174     return QPMS_SUCCESS;
 175   else {
 176     if ( NULL == (c->akn = realloc(c->akn, sizeof(double) * akn_index(lMax + 2, 0))))
 177       abort();
 178     //if ( NULL == (c->bkn = realloc(c->bkn, sizeof(complex double) * bkn_index(lMax + 1, 0))))
 179     //  abort();
 180     for(qpms_l_t n = c->lMax+1; n <= lMax + 1; ++n)
 181       for(qpms_l_t k = 0; k <= n; ++k)
 182         c->akn[akn_index(n,k)] = exp(lgamma(n + k + 1) - k*M_LN2 - lgamma(k + 1) - lgamma(n - k + 1));
 183     // ... TODO derivace
 184     c->lMax = lMax;
 185     return QPMS_SUCCESS;
 186   }
 187 }
 188
 189 complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, complex double x) {
 190   if(QPMS_SUCCESS != qpms_sbessel_calculator_ensure_lMax(c, n))
 191     abort();
 192   complex double z = I/x;
 193   complex double result = 0;
 194   for (qpms_l_t k = n; k >= 0; --k)
 195     // can we use fma for complex?
 196     //result = fma(result, z, c->akn(n, k));
 197     result = result * z + c->akn[akn_index(n,k)];
 198   result *= z * ipow(-n-2) * cexp(I * x);
 199   return result;
 200 }
 201
 202 qpms_errno_t qpms_sbessel_calc_h1_fill(qpms_sbessel_calculator_t * const c,
 203     const qpms_l_t lMax, const complex double x, complex double * const target) {
 204   if(QPMS_SUCCESS != qpms_sbessel_calculator_ensure_lMax(c, lMax))
 205     abort();
 206   memset(target, 0, sizeof(complex double) * lMax);
 207   complex double kahancomp[lMax];
 208   memset(kahancomp, 0, sizeof(complex double) * lMax);
 209   for(qpms_l_t k = 0; k <= lMax; ++k){
 210     double xp = cpow(x, -k-1);
 211     for(qpms_l_t l = k; l <= lMax; ++l)
 212       ckahanadd(target + l, kahancomp + l, c->akn[akn_index(l,k)] * xp * ipow(k-l-1));
 213   }
 214   complex double eix = cexp(I * x);
 215   for (qpms_l_t l = 0; l <= lMax; ++l)
 216     target[l] *= eix;
 217   return QPMS_SUCCESS;
 218 }
 219
 220 qpms_sbessel_calculator_t *qpms_sbessel_calculator_init() {
 221   qpms_sbessel_calculator_t *c = malloc(sizeof(qpms_sbessel_calculator_t));
 222   c->akn = NULL;
 223   //c->bkn = NULL;
 224   c->lMax = -1;
 225   return c;
 226 }
 227
 228 void qpms_sbessel_calculator_pfree(qpms_sbessel_calculator_t *c) {
 229   if(c->akn) free(c->akn);
 230   //if(c->bkn) free(c->bkn);
 231   free(c);
 232 }