oldtests/bessel_precision/besseltest_orig.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <gsl/gsl_sf_bessel.h>
   5
   6 int main() {
   7         int lMax;
   8         while (1 == scanf("%d", &lMax)){
   9                 double x;
  10                 if (1 != scanf(" %lf", &x))
  11                         abort();
  12                 double orig[lMax+1], gsl[lMax+2], relerr[lMax+1];
  13                 for(int l=0; l <= lMax; ++l)
  14                         if (1 != scanf(" %lf", orig+l))
  15                                 abort();
  16 #if defined JTEST || defined DJTEST
  17                 if(gsl_sf_bessel_jl_array(lMax+1,x,gsl))
  18 #elif defined YTEST || defined DYTEST
  19                 if(gsl_sf_bessel_yl_array(lMax+1,x,gsl))
  20 #else
  21                 if(gsl_sf_bessel_jl_steed_array(lMax+1,x,gsl))
  22 #endif
  23                         abort();
  24 #if defined DJTEST || defined DYTEST || defined DJTEST_STEED
  25 #if 1
  26                 for (int l = 0; l <= lMax; ++l)
  27                         gsl[l] = -gsl[l+1] + (l/x) * gsl[l];
  28 #else
  29         // todo varianta 10.51.2
  30 #endif
  31 #endif
  32                 for (int l = 0; l <= lMax; ++l)
  33                         relerr[l] = fabs(gsl[l] - orig[l]) / (fabs(gsl[l]) + fabs(orig[l])) * 2;
  34                 printf("x = %.16g\n", x);
  35                 printf("orig: ");
  36                 for (int l = 0; l <= lMax; ++l)
  37                         printf("%.16g ", orig[l]);
  38                 printf("\ngsl:  ");
  39                 for (int l = 0; l <= lMax; ++l)
  40                         printf("%.16g ", gsl[l]);
  41                 printf("\nrerr:  ");
  42                 for (int l = 0; l <= lMax; ++l)
  43                         printf("%.16g ", relerr[l]);
  44                 putchar('\n');
  45         }
  46 }
  47
  48