Various uncommitted mods and fixes

Former-commit-id: a9c5e0b73d23eae4b8b5fe5223f21a391d321b1e

qpms/bessel.c

@@ -2,8 +2,10 @@
 #include "qpms_specfunc.h"
 #include <stdlib.h>
 #include <stddef.h>
+#include <string.h>
 #include "kahansum.h"
 #include <gsl/gsl_sf_bessel.h>
+#include <complex.h>
 
 #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 = result * z + c->akn[akn_index(n,k)];
-	result *= z * ipow TODOTODOTODO
+	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_pfree(qpms_sbessel_calculator_t *c) {
-
-void qpms_sbessel_calculator_free_p(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);
 }

qpms/qpms_specfunc.h

@@ -12,16 +12,20 @@ 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);
+qpms_sbessel_calculator_t *qpms_sbessel_calculator_init();
-void qpms_bessel_calculator_free(qpms_bessel_calculator_t *c);
+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);
 
 

tests/bessel_precision/besselDJcases.REMOVED.git-id

@@ -1 +1 @@
-12584d0a62798e0d69d43fc030ba6fad056477e4
+2ffa59b967c465cc97b7e163a563e90148b685ef

tests/bessel_precision/besselDYcases.REMOVED.git-id

@@ -1 +1 @@
-12584d0a62798e0d69d43fc030ba6fad056477e4
+ad2823c93ce22441c2a4b9629299c2ed551cec04

tests/bessel_precision/besselJcases.REMOVED.git-id

@@ -1 +1 @@
-a8fbb5a38d9893cec329d03ffe9d8426203b9043
+89b9e7bf3a17acbb8d59501eed51009ff1a8b202

tests/bessel_precision/besselYcases.REMOVED.git-id

@@ -1 +1 @@
-5d2526ef87d9a1b5f5660fea4e9c5ccf1ce9394d
+019a7772dac65ea5e0bc48bfd2abca1deefa9c34

tests/bessel_precision/besseltest.c

@@ -2,28 +2,42 @@
 #include <stdio.h>
 #include <math.h>
 
+#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 defined JTEST_QPMS || defined DJTEST_QPMS || defined YTEST_QPMS || defined DYTEST_QPMS
-		if (gsl_sf_bessel_jl_array(lMax+1, x, gsl))
+		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))
+		if (gsl_sf_bessel_yl_array(lMax+1, x, gsl)) abort();
-#else
+#elif defined JTEST_STEED || DJTEST_STEED
-		if (gsl_sf_bessel_jl_steed_array(lMax+1, x, gsl))
+		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;
 }
 

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
+
 

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))