[gnumeric] BesselI: Use series for small x.

[Morten Welinder <mortenw src gnome org>]

commit c57ea4f52355dabab18e35d2204dce4f25d2533f
Author: Morten Welinder <terra gnome org>
Date:   Fri Jan 29 21:19:41 2016 -0500

    BesselI: Use series for small x.

 ChangeLog       |    4 ++++
 src/sf-bessel.c |   25 ++++++++++++++++---------
 2 files changed, 20 insertions(+), 9 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 14a1d55..8223092 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,9 @@
 2016-01-29  Morten Welinder  <terra gnome org>
 
+       * src/sf-bessel.c (gnm_bessel_i): Use series for small x.
+       (bessel_ij_series): Rename from bessel_j_series and extend to
+       handle i also.
+
        * src/sf-gamma.c (qbetaf): Improve accuracy in the case where one
        argument is less than one.
 
diff --git a/src/sf-bessel.c b/src/sf-bessel.c
index 2861c2a..e6c106e 100644
--- a/src/sf-bessel.c
+++ b/src/sf-bessel.c
@@ -1158,7 +1158,7 @@ L420:
 /* ------------------------------------------------------------------------ */
 
 static gboolean
-bessel_j_series_domain (gnm_float x, gnm_float v)
+bessel_ij_series_domain (gnm_float x, gnm_float v)
 {
        // The taylor series is valid for all possible values of x and v,
        // but it isn't efficient and practical for all.
@@ -1185,7 +1185,7 @@ bessel_j_series_domain (gnm_float x, gnm_float v)
 
 
 static gnm_float
-bessel_j_series (gnm_float x, gnm_float v)
+bessel_ij_series (gnm_float x, gnm_float v, gboolean qj)
 {
        GnmQuad qv, qxh, qfv, qs, qt;
        int efv;
@@ -1232,8 +1232,9 @@ bessel_j_series (gnm_float x, gnm_float v)
                        gnm_float t;
 
                        gnm_quad_mul (&qt, &qt, &qxh2);
-                       gnm_quad_init (&qa, -k);
-                       gnm_quad_sub (&qb, &qv, &qa);
+                       gnm_quad_init (&qa, k);
+                       gnm_quad_add (&qb, &qv, &qa);
+                       gnm_quad_init (&qa, qj ? -k : k);
                        gnm_quad_mul (&qa, &qa, &qb);
                        gnm_quad_div (&qt, &qt, &qa);
                        t = gnm_quad_value (&qt);
@@ -2263,8 +2264,8 @@ 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_j_series (x, nu);
-       gnm_float Jmnu = bessel_j_series (x, -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);
        return Ynu;
@@ -2291,11 +2292,11 @@ static void
 hankel1_A1 (gnm_complex *res, gnm_float x, gnm_float nu)
 {
        gnm_float rnu = gnm_floor (nu + 0.5);
-       gnm_float Jnu = bessel_j_series (x, nu);
+       gnm_float Jnu = bessel_ij_series (x, nu, TRUE);
        gnm_float Ynu;
 
        if (gnm_abs (nu - rnu) > 5e-4) {
-               gnm_float Jmnu = bessel_j_series (x, -nu);
+               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;
@@ -2405,7 +2406,7 @@ hankel1 (gnm_complex *res, gnm_float x, gnm_float nu)
                // Algorithm A
                // Deviation: we use the series on a wider domain as our
                // series code uses quad precision.
-               if (bessel_j_series_domain (x, nu))
+               if (bessel_ij_series_domain (x, nu))
                        hankel1_A1 (res, x, nu);
                else if (nu > x && g > 1.5)
                        hankel1_A2 (res, x, nu);
@@ -2421,6 +2422,12 @@ hankel1 (gnm_complex *res, gnm_float x, gnm_float nu)
 gnm_float
 gnm_bessel_i (gnm_float x, gnm_float alpha)
 {
+       if (gnm_isnan (x) || gnm_isnan (alpha))
+               return x + alpha;
+
+       if (bessel_ij_series_domain (x, alpha))
+               return bessel_ij_series (x, alpha, FALSE);
+
        if (x < 0) {
                if (alpha != gnm_floor (alpha))
                        return gnm_nan;