ENH: add changes from Amos (1995) to amos.h

[JaRoSchm]

Reference issue

Closes #137

What does this implement/fix?

This adds the improvements from Amos' 1995 paper to the Amos library. The original Fortran code, which I translated by hand, is available here (functional changes are limited to besy, the old version has been renamed to besyh equivalent to the Fortran code). I became aware of this as I stumbled across this repo, which compares the accuracy of the original Fortran code, scipy, and a rust translation. The deviations for the scaled function yve are fixed by this PR (I suspect the remaining ones to appear due to translation errors from the original translation). Here I reproduced their plot with this PR (the version of scipy in the legend is wrong):

Locally, the xsref tests for yv and jv (which uses yv here) fail on this PR. I think this can be fixed by slightly increasing the respective relative tolerances. I computed the percentiles for the test cases from xsref for yv, yve, jv, and jve for real and complex inputs to show that there are no significant negative impacts of this PR. However, there are not many positive changes either as the regions of improved accuracy seem to not be included in xsref:

yv dd->d (7780 inputs)

	Min	10%	20%	30%	40%	50%	60%	70%	80%	90%	Max
new	0.00e+00	0.00e+00	0.00e+00	0.00e+00	0.00e+00	1.67e-16	3.32e-16	6.12e-16	1.39e-15	1.78e-14	inf
old	0.00e+00	0.00e+00	0.00e+00	0.00e+00	0.00e+00	1.51e-16	2.54e-16	4.57e-16	9.75e-16	8.59e-15	inf

Number of increased relative errors: 1107
Number of decreased relative errors: 505

yv dD->D (17158 inputs)

	Min	10%	20%	30%	40%	50%	60%	70%	80%	90%	Max
new	0.00e+00	0.00e+00	1.33e-16	2.64e-16	5.55e-16	1.59e-15	6.29e-15	1.69e-14	nan	nan	inf
old	0.00e+00	0.00e+00	8.20e-17	2.10e-16	4.34e-16	1.22e-15	5.64e-15	1.54e-14	nan	nan	inf

Number of increased relative errors: 5154
Number of decreased relative errors: 2549

yve dd->d (20 inputs)

	Min	10%	20%	30%	40%	50%	60%	70%	80%	90%	Max
new	0.00e+00	0.00e+00	0.00e+00	0.00e+00	0.00e+00	5.57e-17	1.16e-16	1.45e-16	2.14e-16	3.33e-16	7.71e-16
old	0.00e+00	0.00e+00	0.00e+00	0.00e+00	0.00e+00	0.00e+00	1.23e-16	1.53e-16	2.14e-16	2.87e-16	3.09e-16

Number of increased relative errors: 4
Number of decreased relative errors: 1

yve dD->D (66 inputs)

	Min	10%	20%	30%	40%	50%	60%	70%	80%	90%	Max
new	0.00e+00	5.35e-17	8.56e-17	1.07e-16	1.15e-16	1.67e-16	2.98e-16	4.07e-16	4.20e-16	4.25e-16	1.29e-15
old	1.00e-16	1.15e-16	2.56e-16	2.85e-16	3.30e-16	3.52e-16	3.55e-16	4.25e-16	5.09e-16	8.56e-16	1.29e-15

Number of increased relative errors: 13
Number of decreased relative errors: 44

jv dd->d (3416 inputs)

	Min	10%	20%	30%	40%	50%	60%	70%	80%	90%	Max
new	0.00e+00	0.00e+00	1.42e-16	3.16e-16	5.95e-16	1.33e-15	2.42e-15	4.89e-15	9.30e-15	2.14e-14	inf
old	0.00e+00	0.00e+00	1.40e-16	3.16e-16	5.95e-16	1.33e-15	2.42e-15	4.91e-15	9.28e-15	2.14e-14	inf

Number of increased relative errors: 56
Number of decreased relative errors: 57

jv dD->D (3374 inputs)

	Min	10%	20%	30%	40%	50%	60%	70%	80%	90%	Max
new	0.00e+00	1.62e-16	3.18e-16	6.69e-16	1.15e-15	2.21e-15	3.56e-15	6.45e-15	1.15e-14	2.46e-14	inf
old	0.00e+00	1.59e-16	3.16e-16	6.48e-16	1.12e-15	2.14e-15	3.54e-15	6.33e-15	1.12e-14	2.44e-14	inf

Number of increased relative errors: 215
Number of decreased relative errors: 69

jve dd->d (21 inputs)

	Min	10%	20%	30%	40%	50%	60%	70%	80%	90%	Max
new	0.00e+00	1.26e-16	1.26e-16	1.34e-16	1.34e-16	1.97e-16	1.97e-16	2.15e-16	2.79e-16	3.67e-15	3.67e-15
old	0.00e+00	1.26e-16	1.26e-16	1.34e-16	1.34e-16	1.97e-16	1.97e-16	2.15e-16	2.79e-16	3.67e-15	3.67e-15

Number of increased relative errors: 0
Number of decreased relative errors: 0

jve dD->D (66 inputs)

	Min	10%	20%	30%	40%	50%	60%	70%	80%	90%	Max
new	0.00e+00	1.15e-16	1.27e-16	1.74e-16	2.36e-16	2.67e-16	2.98e-16	4.25e-16	4.42e-16	1.23e-15	1.79e-15
old	0.00e+00	1.15e-16	1.27e-16	1.74e-16	2.36e-16	2.67e-16	2.98e-16	4.25e-16	4.42e-16	1.23e-15	1.79e-15

Number of increased relative errors: 0
Number of decreased relative errors: 0

Inputs significantly improved by this PR:
$\nu$ = 100: z = +-100 +- 0.001/0.1/1/5/10/30/50i, z = +- 200/300 +- 0.001/0.1/1/5/10/30/50/100i, ...
$\nu$ = 500, z = +-1000 +- 0.001/0.1/1/5/10/30/50/100/200/500i, ...

Additional information

None

AI Generation Disclosure

No AI tools used

[JaRoSchm]

Actually, all the other outliers for kve, hankel1e, and hankel2e are fixed by some combination of #92, #93, and #132. I did not apply them one by one, but if all of them are applied we have:

[JaRoSchm]

Maybe, we should use the new implementation for yve with complex z only and continue to use the old one for the other seven cases?