oldtests/bessel_precision/besseltest.c

   1 #include <gsl/gsl_sf_bessel.h>
   2 #include <stdio.h>
   3 #include <math.h>
   4
   5 #if defined JTEST_QPMS || defined DJTEST_QPMS || defined YTEST_QPMS || defined DYTEST_QPMS
   6 #include "../../qpms/qpms_specfunc.h"
   7 #endif
   8
   9
  10 int main() {
  11         int lMax;
  12 #if defined JTEST_QPMS || defined DJTEST_QPMS || defined YTEST_QPMS || defined DYTEST_QPMS
  13         qpms_sbessel_calculator_t *c = qpms_sbessel_calculator_init();
  14 #endif
  15         while (1 == scanf("%d", &lMax)) {
  16                 double x;
  17                 if (1 != scanf("%lf", &x))
  18                         abort();
  19                 double gsl[lMax+2], relerr[lMax+1], orig[lMax+1];
  20
  21                 for (int l = 0; l <= lMax; l++)
  22                         if (1 != scanf("%lf", orig+l))
  23                                 abort();
  24 #if defined JTEST_QPMS || defined DJTEST_QPMS || defined YTEST_QPMS || defined DYTEST_QPMS
  25                 complex double hankel[lMax+2];
  26                 qpms_sbessel_calc_h1_fill(c, lMax+1, x, hankel);
  27                 for(int l = 0; l <= lMax+1; l++)
  28         #if defined JTEST_QPMS
  29                         gsl[l] = creal(hankel[l]);
  30         #endif
  31 #elif defined JTEST || defined DJTEST
  32                 if (gsl_sf_bessel_jl_array(lMax+1, x, gsl)) abort();
  33 #elif defined YTEST || defined DYTEST
  34                 if (gsl_sf_bessel_yl_array(lMax+1, x, gsl)) abort();
  35 #elif defined JTEST_STEED || DJTEST_STEED
  36                 if (gsl_sf_bessel_jl_steed_array(lMax+1, x, gsl)) abort();
  37 #endif
  38 #if defined DJTEST || defined DYTEST || defined DJTEST_STEED
  39                 for (int l = 0; l <= lMax; l++)
  40                         gsl[l] = -gsl[l+1] + ((double)l/x) * gsl[l];
  41 #endif
  42                 printf("x = %.16g, lMax = %d:\nsage: ", x, lMax);
  43                 for (int l = 0; l <= lMax; l++)
  44                         printf("%.16g ", orig[l]);
  45                 printf("\ngsl:  ");
  46                 for (int l = 0; l <= lMax; l++)
  47                         printf("%.16g ", gsl[l]);
  48                 printf("\nrell: ");
  49                 for (int l = 0; l <= lMax; l++)
  50                         printf("%.16g ", 2 * fabs(gsl[l] - orig[l]) / (fabs(gsl[l]) + fabs(gsl[l])));
  51                 putchar('\n');
  52         }
  53 #if defined JTEST_QPMS || defined DJTEST_QPMS || defined YTEST_QPMS || defined DYTEST_QPMS
  54         qpms_sbessel_calculator_pfree(c);
  55 #endif
  56         return 0;
  57 }
  58