diff --git a/qpms/bessel.c b/qpms/bessel.c index 0cbc45d..78097a7 100644 --- a/qpms/bessel.c +++ b/qpms/bessel.c @@ -2,8 +2,10 @@ #include "qpms_specfunc.h" #include #include +#include #include "kahansum.h" #include +#include #ifndef M_LN2 #define M_LN2 0.69314718055994530942 /* log_e 2 */ @@ -64,46 +66,56 @@ static inline qpms_errno_t qpms_sbessel_calculator_ensure_lMax(qpms_sbessel_calc abort(); //if ( NULL == (c->bkn = realloc(c->bkn, sizeof(complex double) * bkn_index(lMax + 1, 0)))) // abort(); - for(qpms_l_t n = c->lMax+1; l <= lMax + 1; ++l) + for(qpms_l_t n = c->lMax+1; n <= lMax + 1; ++n) for(qpms_l_t k = 0; k <= n; ++k) c->akn[akn_index(n,k)] = exp(lgamma(n + k + 1) - k*M_LN2 - lgamma(k + 1) - lgamma(n - k + 1)); // ... TODO derivace c->lMax = lMax; - return QPMS_SUCCESS + return QPMS_SUCCESS; } } complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, double x) { - double imaginary z = I/x; + if(QPMS_SUCCESS != qpms_sbessel_calculator_ensure_lMax(c, n)) + abort(); + complex double z = I/x; // FIXME this should be imaginary double, but gcc is broken? complex double result = 0; for (qpms_l_t k = n; k >= 0; --k) // can we use fma for complex? //result = fma(result, z, c->akn(n, k)); - result = result * z + c->akn(n,k); - result *= z * ipow TODOTODOTODO - - - - - - - - - + result = result * z + c->akn[akn_index(n,k)]; + result *= z * ipow(-n-2) * cexp(I * x); + return result; +} +qpms_errno_t qpms_sbessel_calc_h1_fill(qpms_sbessel_calculator_t * const c, + const qpms_l_t lMax, const double x, complex double * const target) { + if(QPMS_SUCCESS != qpms_sbessel_calculator_ensure_lMax(c, lMax)) + abort(); + memset(target, 0, sizeof(complex double) * lMax); + complex double kahancomp[lMax]; + memset(kahancomp, 0, sizeof(complex double) * lMax); + for(qpms_l_t k = 0; k <= lMax; ++k){ + double xp = pow(x, -k-1); + for(qpms_l_t l = k; l <= lMax; ++l) + ckahanadd(target + l, kahancomp + l, c->akn[akn_index(l,k)] * xp * ipow(k-l-1)); + } + complex double eix = cexp(I * x); + for (qpms_l_t l = 0; l <= lMax; ++l) + target[l] *= eix; + return QPMS_SUCCESS; +} qpms_sbessel_calculator_t *qpms_sbessel_calculator_init() { qpms_sbessel_calculator_t *c = malloc(sizeof(qpms_sbessel_calculator_t)); c->akn = NULL; - c->bkn = NULL; + //c->bkn = NULL; c->lMax = -1; return c; } - - -void qpms_sbessel_calculator_free_p(qpms_sbessel_calculator_t *c) { +void qpms_sbessel_calculator_pfree(qpms_sbessel_calculator_t *c) { if(c->akn) free(c->akn); - if(c->bkn) free(c->bkn); + //if(c->bkn) free(c->bkn); free(c); } diff --git a/qpms/qpms_specfunc.h b/qpms/qpms_specfunc.h index 064d46b..8771339 100644 --- a/qpms/qpms_specfunc.h +++ b/qpms/qpms_specfunc.h @@ -12,18 +12,22 @@ qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, double x, co -typedef struct qpms_bessel_calculator_t { +typedef struct { qpms_l_t lMax; double *akn; // coefficients as in DLMF 10.49.1 //complex double *bkn; // coefficients of the derivatives -}; +} qpms_sbessel_calculator_t; -qpms_bessel_calculator_t *qpms_bessel_calculator_init(qpms_l_t lMax); -void qpms_bessel_calculator_free(qpms_bessel_calculator_t *c); +qpms_sbessel_calculator_t *qpms_sbessel_calculator_init(); +void qpms_sbessel_calculator_pfree(qpms_sbessel_calculator_t *c); -qpms_errno_t qpms_sph_bessel_calc_fill(qpms_bessel_calculator_t *c, qpms_bessel_t typ, qpms_l_t lmax, +qpms_errno_t qpms_sbessel_calc_fill(qpms_sbessel_calculator_t *c, qpms_bessel_t typ, qpms_l_t lmax, double x, complex double *result_array); +complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, double x); +qpms_errno_t qpms_sbessel_calc_h1_fill(qpms_sbessel_calculator_t *c, qpms_l_t lmax, + double x, complex double *result_array); + /****************************************************************************** * Legendre functions and their "angular derivatives" * diff --git a/tests/bessel_precision/besselDJcases.REMOVED.git-id b/tests/bessel_precision/besselDJcases.REMOVED.git-id index 05d7222..d6aa9cf 100644 --- a/tests/bessel_precision/besselDJcases.REMOVED.git-id +++ b/tests/bessel_precision/besselDJcases.REMOVED.git-id @@ -1 +1 @@ -12584d0a62798e0d69d43fc030ba6fad056477e4 \ No newline at end of file +2ffa59b967c465cc97b7e163a563e90148b685ef \ No newline at end of file diff --git a/tests/bessel_precision/besselDYcases.REMOVED.git-id b/tests/bessel_precision/besselDYcases.REMOVED.git-id index 05d7222..69401ea 100644 --- a/tests/bessel_precision/besselDYcases.REMOVED.git-id +++ b/tests/bessel_precision/besselDYcases.REMOVED.git-id @@ -1 +1 @@ -12584d0a62798e0d69d43fc030ba6fad056477e4 \ No newline at end of file +ad2823c93ce22441c2a4b9629299c2ed551cec04 \ No newline at end of file diff --git a/tests/bessel_precision/besselJcases.REMOVED.git-id b/tests/bessel_precision/besselJcases.REMOVED.git-id index 9b87a7e..a81e813 100644 --- a/tests/bessel_precision/besselJcases.REMOVED.git-id +++ b/tests/bessel_precision/besselJcases.REMOVED.git-id @@ -1 +1 @@ -a8fbb5a38d9893cec329d03ffe9d8426203b9043 \ No newline at end of file +89b9e7bf3a17acbb8d59501eed51009ff1a8b202 \ No newline at end of file diff --git a/tests/bessel_precision/besselYcases.REMOVED.git-id b/tests/bessel_precision/besselYcases.REMOVED.git-id index 03b21e5..75ce606 100644 --- a/tests/bessel_precision/besselYcases.REMOVED.git-id +++ b/tests/bessel_precision/besselYcases.REMOVED.git-id @@ -1 +1 @@ -5d2526ef87d9a1b5f5660fea4e9c5ccf1ce9394d \ No newline at end of file +019a7772dac65ea5e0bc48bfd2abca1deefa9c34 \ No newline at end of file diff --git a/tests/bessel_precision/besseltest.c b/tests/bessel_precision/besseltest.c index 794456a..a4f93bf 100644 --- a/tests/bessel_precision/besseltest.c +++ b/tests/bessel_precision/besseltest.c @@ -2,28 +2,42 @@ #include #include +#if defined JTEST_QPMS || defined DJTEST_QPMS || defined YTEST_QPMS || defined DYTEST_QPMS +#include "../../qpms/qpms_specfunc.h" +#endif + int main() { int lMax; +#if defined JTEST_QPMS || defined DJTEST_QPMS || defined YTEST_QPMS || defined DYTEST_QPMS + qpms_sbessel_calculator_t *c = qpms_sbessel_calculator_init(); +#endif while (1 == scanf("%d", &lMax)) { double x; if (1 != scanf("%lf", &x)) abort(); double gsl[lMax+2], relerr[lMax+1], orig[lMax+1]; + for (int l = 0; l <= lMax; l++) if (1 != scanf("%lf", orig+l)) abort(); -#if defined JTEST || defined DJTEST - if (gsl_sf_bessel_jl_array(lMax+1, x, gsl)) +#if defined JTEST_QPMS || defined DJTEST_QPMS || defined YTEST_QPMS || defined DYTEST_QPMS + complex double hankel[lMax+2]; + qpms_sbessel_calc_h1_fill(c, lMax+1, x, hankel); + for(int l = 0; l <= lMax+1; l++) + #if defined JTEST_QPMS + gsl[l] = creal(hankel[l]); + #endif +#elif defined JTEST || defined DJTEST + if (gsl_sf_bessel_jl_array(lMax+1, x, gsl)) abort(); #elif defined YTEST || defined DYTEST - if (gsl_sf_bessel_yl_array(lMax+1, x, gsl)) -#else - if (gsl_sf_bessel_jl_steed_array(lMax+1, x, gsl)) + if (gsl_sf_bessel_yl_array(lMax+1, x, gsl)) abort(); +#elif defined JTEST_STEED || DJTEST_STEED + if (gsl_sf_bessel_jl_steed_array(lMax+1, x, gsl)) abort(); #endif - abort(); #if defined DJTEST || defined DYTEST || defined DJTEST_STEED for (int l = 0; l <= lMax; l++) - gsl[l] = -gsl[l+1] + (l/x) * gsl[l]; + gsl[l] = -gsl[l+1] + ((double)l/x) * gsl[l]; #endif printf("x = %.16g, lMax = %d:\nsage: ", x, lMax); for (int l = 0; l <= lMax; l++) @@ -36,6 +50,9 @@ int main() { printf("%.16g ", 2 * fabs(gsl[l] - orig[l]) / (fabs(gsl[l]) + fabs(gsl[l]))); putchar('\n'); } +#if defined JTEST_QPMS || defined DJTEST_QPMS || defined YTEST_QPMS || defined DYTEST_QPMS + qpms_sbessel_calculator_pfree(c); +#endif return 0; } diff --git a/tests/bessel_precision/compileall.sh b/tests/bessel_precision/compileall.sh index bcac297..1d7b2c5 100755 --- a/tests/bessel_precision/compileall.sh +++ b/tests/bessel_precision/compileall.sh @@ -6,4 +6,6 @@ c99 -ggdb -o jtest_steed -DJTEST_STEED besseltest.c -lgsl -lblas -lm c99 -ggdb -o djtest -DDJTEST besseltest.c -lgsl -lblas -lm c99 -ggdb -o dytest -DDYTEST besseltest.c -lgsl -lblas -lm c99 -ggdb -o djtest_steed -DDJTEST_STEED besseltest.c -lgsl -lblas -lm +c99 -ggdb -o jtest_qpms -DJTEST_QPMS besseltest.c ../../qpms/bessel.c -lgsl -lm -lblas + diff --git a/tests/bessel_precision/sagebesselgen.py b/tests/bessel_precision/sagebesselgen.py index d232f0d..12a6e91 100644 --- a/tests/bessel_precision/sagebesselgen.py +++ b/tests/bessel_precision/sagebesselgen.py @@ -20,6 +20,9 @@ def printbesselYrow(lMax, x, file=sys.stdout): print(N(spherical_bessel_Y(l,x)), end = ' ', file=file) print('', file=file) + +#cf DLMF 10.51.2 + def printbesselDJrow(lMax, x, file=sys.stdout): print(lMax, N(x), end=' ', file=file); for l in range(lMax+1): @@ -36,6 +39,9 @@ def printbesselDYrow(lMax, x, file=sys.stdout): print('', file=file) +def ank(n, k): + return factorial(n+k)/2**k/factorial(k)/factorial(n-k) + def genall(lMax): f = open('besselDJcases', 'w') for o in IntegerRange(1,100): @@ -58,3 +64,16 @@ def genall(lMax): printbesselYrow(lMax, 1/o, file=f) printbesselYrow(lMax, o/sqrt(3), file=f) + +import math +M_LN2 = math.log(2) + +def ankf(n,k): + n = float(n) + k = float(k) + return math.exp(math.lgamma(n+k+1) - k * M_LN2 - math.lgamma(k+1) - math.lgamma(n-k+1)) + +def ankrelerr(n,k): + a = ank(n,k) + b = ankf(n,k) + return 2 * abs(a - b)/(abs(a)+abs(b))