[gnumeric] BESSELJ: Use taylor series when appropriate.

[Morten Welinder <mortenw src gnome org>]

commit f2966e89259db55b9f039602cdded11bac24cd7f
Author: Morten Welinder <terra gnome org>
Date:   Sun Dec 1 18:43:27 2013 -0500

    BESSELJ: Use taylor series when appropriate.

 ChangeLog       |    3 +++
 NEWS            |    1 +
 src/sf-bessel.c |   26 ++++++++++++++++++++++++++
 3 files changed, 30 insertions(+), 0 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index bb0ced7..e89879d 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,8 @@
 2013-12-01  welinder  <terra gnome org>
 
+       * src/sf-bessel.c (bessel_j): Use the taylor series in the
+       parameter range where that makes sense.
+
        * src/sf-gamma.c (gnm_lbeta3): Improve accuracy.
        (gnm_beta): Ditto.
 
diff --git a/NEWS b/NEWS
index 12124f0..c3a3f87 100644
--- a/NEWS
+++ b/NEWS
@@ -6,6 +6,7 @@ Andreas:
 Morten:
        * Extend POCHHAMMER to negative values.
        * Improve accuracy of BETA and BETALN.
+       * Improve accuracy of BESSELJ.
 
 --------------------------------------------------------------------------
 Gnumeric 1.12.9
diff --git a/src/sf-bessel.c b/src/sf-bessel.c
index c3233bc..97cb122 100644
--- a/src/sf-bessel.c
+++ b/src/sf-bessel.c
@@ -11,6 +11,8 @@
 #define MATHLIB_WARNING2 g_warning
 #define MATHLIB_WARNING4 g_warning
 
+static gnm_float bessel_j_series (gnm_float x, gnm_float alpha);
+
 static gnm_float bessel_y_ex(gnm_float x, gnm_float alpha, gnm_float *by);
 static gnm_float bessel_k(gnm_float x, gnm_float alpha, gnm_float expo);
 static gnm_float bessel_y(gnm_float x, gnm_float alpha);
@@ -683,6 +685,9 @@ static gnm_float bessel_j(gnm_float x, gnm_float alpha)
               ((alpha == na) ? 0 :
               bessel_y(x, -alpha) * gnm_sinpi(alpha)));
     }
+    if (x < 2 || x * x / 4 < alpha)
+           return bessel_j_series (x, alpha);
+
     nb = 1 + (long)na; /* nb-1 <= alpha < nb */
     alpha -= (gnm_float)(nb-1);
 #ifdef MATHLIB_STANDALONE
@@ -2230,6 +2235,27 @@ L450:
 
 /* ------------------------------------------------------------------------ */
 
+static gnm_float
+bessel_j_series (gnm_float x, gnm_float alpha)
+{
+       gnm_float s, t;
+       gnm_float xh = x / 2, xh2 = xh * xh;
+       int k;
+
+       t = gnm_pow (xh, alpha) / gnm_fact (alpha);
+       s = t;
+       for (k = 1; k < 1000; k++) {
+               t *= -xh2 / (k * (k + alpha));
+               s += t;
+               if (k > 5 && gnm_abs (t) < GNM_EPSILON * gnm_abs (s))
+                       break;
+       }
+
+       return s;
+}
+
+/* ------------------------------------------------------------------------ */
+
 gnm_float
 gnm_bessel_i (gnm_float x, gnm_float alpha)
 {