From fb9bbff4b8abb95db1042d709b06cba99e5007a4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Wed, 14 Mar 2018 14:10:11 +0200
Subject: [PATCH] gsl_sf_bessel_jl_array and gsl_sf_bessel_jl_steed_array are
 inprecise.

Must do a different implementation.


Former-commit-id: 37da4cf29b5d8706c639a9a3af61653d90ccc37e
---
 notes/tftests.lyx                    | 17 ++++++++++--
 tests/bessel_precision/besseltest.c  | 41 ++++++++++++++++++++++++++++
 tests/bessel_precision/compileall.sh |  3 +-
 3 files changed, 57 insertions(+), 4 deletions(-)
 create mode 100644 tests/bessel_precision/besseltest.c

diff --git a/notes/tftests.lyx b/notes/tftests.lyx
index 3eafdef..6adcfbc 100644
--- a/notes/tftests.lyx
+++ b/notes/tftests.lyx
@@ -343,7 +343,7 @@ key "xu_efficient_1998"
 \begin_inset Formula $B_{mn\mu\nu}$
 \end_inset
 
- se liší v několika ohledech:
+ se lišejí v několika ohledech:
 \end_layout
 
 \begin_layout Enumerate
@@ -389,9 +389,18 @@ Ve starší práci je poslední člen sumy
 \begin_inset Formula $\left(m,n,\mu,\nu\right)=\left(0,1,-1,1\right)$
 \end_inset
 
-!!! (Opraveno, použit Cruzanův tvar v 
+!!! (Při numerickém srovnání Xuových vzorců 
 \begin_inset CommandInset citation
 LatexCommand cite
+after "(77–80)"
+key "xu_calculation_1996"
+
+\end_inset
+
+ ve staré práci a cruzanovských vzorců 
+\begin_inset CommandInset citation
+LatexCommand cite
+after "(59–61, ...)"
 key "xu_efficient_1998"
 
 \end_inset
@@ -430,7 +439,9 @@ Legendreovy polynomy,
 \end_layout
 
 \begin_layout Itemize
-Besselovy funkce.
+Besselovy funkce – nepřesné jak sviňa zejména u derivací besselových funkcí
+ prvního druhu.
+ Nutno zvolit jinou implementaci.
 \end_layout
 
 \end_deeper
diff --git a/tests/bessel_precision/besseltest.c b/tests/bessel_precision/besseltest.c
new file mode 100644
index 0000000..794456a
--- /dev/null
+++ b/tests/bessel_precision/besseltest.c
@@ -0,0 +1,41 @@
+#include <gsl/gsl_sf_bessel.h>
+#include <stdio.h>
+#include <math.h>
+
+
+int main() {
+	int lMax;
+	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))
+#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))
+#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];
+#endif
+		printf("x = %.16g, lMax = %d:\nsage: ", x, lMax);
+		for (int l = 0; l <= lMax; l++)
+			printf("%.16g ", orig[l]);
+		printf("\ngsl:  ");
+		for (int l = 0; l <= lMax; l++)
+			printf("%.16g ", gsl[l]);
+		printf("\nrell: ");
+		for (int l = 0; l <= lMax; l++)
+			printf("%.16g ", 2 * fabs(gsl[l] - orig[l]) / (fabs(gsl[l]) + fabs(gsl[l])));
+		putchar('\n');
+	}
+	return 0;
+}
+
diff --git a/tests/bessel_precision/compileall.sh b/tests/bessel_precision/compileall.sh
index e16e440..bcac297 100755
--- a/tests/bessel_precision/compileall.sh
+++ b/tests/bessel_precision/compileall.sh
@@ -2,7 +2,8 @@
 
 c99 -ggdb -o jtest -DJTEST besseltest.c -lgsl -lblas -lm
 c99 -ggdb -o ytest -DYTEST besseltest.c -lgsl -lblas -lm
+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 -DJTEST_STEED besseltest.c -lgsl -lblas -lm
+c99 -ggdb -o djtest_steed -DDJTEST_STEED besseltest.c -lgsl -lblas -lm