[gnumeric] BESSELY: Minor improvements for integer order.

[Morten Welinder <mortenw src gnome org>]

commit 61fe7c4f27f517fccaa007a7ab491563c6a5a854
Author: Morten Welinder <terra gnome org>
Date:   Tue Mar 28 14:27:00 2017 -0400

    BESSELY: Minor improvements for integer order.

 ChangeLog       |    6 ++++++
 NEWS            |    1 +
 src/numbers.h   |    2 ++
 src/sf-bessel.c |   30 ++++++++++++++++++++++--------
 4 files changed, 31 insertions(+), 8 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 26a6160..5b2a367 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,9 @@
+2017-03-28  Morten Welinder  <terra gnome org>
+
+       * src/sf-bessel.c (y_via_j_series): Use gnm_yn if we can.
+       (hankel1_A1): Use gnm_yn if we can.  Near integer order, center
+       around the integer with odd count when we can.
+
 2017-03-21  Morten Welinder  <terra gnome org>
 
        * src/sstest.c (rand_fractile_test): Allow 4*sqrt(expected)
diff --git a/NEWS b/NEWS
index e9ef541..9ca860f 100644
--- a/NEWS
+++ b/NEWS
@@ -2,6 +2,7 @@ Gnumeric 1.12.35
 
 Morten:
        * Test suite improvements.
+       * Minor BESSELY improvements.
 
 --------------------------------------------------------------------------
 Gnumeric 1.12.34
diff --git a/src/numbers.h b/src/numbers.h
index e38f632..7d29c52 100644
--- a/src/numbers.h
+++ b/src/numbers.h
@@ -91,6 +91,7 @@ typedef long double gnm_float;
 #define gnm_cos cosl
 #define gnm_tan tanl
 #endif
+#define gnm_yn ynl
 
 #define GNM_FORMAT_e   "Le"
 #define GNM_FORMAT_E   "LE"
@@ -206,6 +207,7 @@ typedef double gnm_float;
 #define gnm_cos cos
 #define gnm_tan tan
 #endif
+#define gnm_yn yn
 
 #define GNM_FORMAT_e   "e"
 #define GNM_FORMAT_E   "E"
diff --git a/src/sf-bessel.c b/src/sf-bessel.c
index 5a6afce..4a6f487 100644
--- a/src/sf-bessel.c
+++ b/src/sf-bessel.c
@@ -2212,12 +2212,20 @@ static gnm_float
 y_via_j_series (gnm_float nu, const gnm_float *args)
 {
        gnm_float x = args[0];
-       gnm_float c = gnm_cospi (nu);
-       gnm_float s = gnm_sinpi (nu);
-       gnm_float c_Jnu = c == 0 ? 0 : c * bessel_ij_series (x, nu, TRUE);
-       gnm_float Jmnu = bessel_ij_series (x, -nu, TRUE);
-       gnm_float Ynu = (c_Jnu - Jmnu) / s;
-       if (0) g_printerr ("via: %.20g %.20g %.20g %.20g %.20g\n", x, nu, c_Jnu, Jmnu, Ynu);
+       gnm_float Ynu;
+       if (nu == gnm_floor (nu) && gnm_abs (nu) < G_MAXINT) {
+               Ynu = gnm_yn ((int)nu, x);
+               if (0) g_printerr ("via: %.20g %.20g %.20g\n", x, nu, Ynu);
+       } else {
+               gnm_float c = gnm_cospi (nu);
+               gnm_float s = gnm_sinpi (nu);
+               gnm_float c_Jnu = c == 0
+                       ? 0
+                       : c * bessel_ij_series (x, nu, TRUE);
+               gnm_float Jmnu = bessel_ij_series (x, -nu, TRUE);
+               Ynu = (c_Jnu - Jmnu) / s;
+               if (0) g_printerr ("via: %.20g %.20g %.20g %.20g %.20g\n", x, nu, c_Jnu, Jmnu, Ynu);
+       }
        return Ynu;
 }
 
@@ -2241,20 +2249,26 @@ hankel1_B2 (gnm_float x, gnm_float nu, size_t N)
 static gnm_complex
 hankel1_A1 (gnm_float x, gnm_float nu)
 {
-       gnm_float rnu = gnm_floor (nu + 0.5);
+       gnm_float rnu = gnm_floor (nu + 0.49); // Close enough
        gnm_float Jnu = bessel_ij_series (x, nu, TRUE);
        gnm_float Ynu;
+       gboolean use_yn = (gnm_abs (rnu) < G_MAXINT - 1);
 
        if (gnm_abs (nu - rnu) > 5e-4) {
                gnm_float Jmnu = bessel_ij_series (x, -nu, TRUE);
                gnm_float c = gnm_cospi (nu);
                gnm_float s = gnm_sinpi (nu);
                Ynu = ((c == 0 ? 0 : Jnu * c) - Jmnu) / s;
+       } else if (use_yn && nu == rnu) {
+               Ynu = gnm_yn ((int)rnu, x);
        } else {
                size_t N = 6;
                gnm_float dnu = 1e-3;
                gnm_float args[1] = { x };
-               Ynu = chebyshev_interpolant (N, rnu - dnu, rnu + dnu, nu,
+               gnm_float nul = rnu - dnu, nur = rnu + dnu;
+               if (use_yn)
+                       N |= 1; // Odd, so we use rnu
+               Ynu = chebyshev_interpolant (N, nul, nur, nu,
                                             y_via_j_series, args);
        }