[gnumeric] BESSEL[IJKY]: Improve accuracy.

[Morten Welinder <mortenw src gnome org>]

commit 4bfc616bf0bd26c9a2a57e8cc817cd207d1f1242
Author: Morten Welinder <terra gnome org>
Date:   Sun Nov 3 15:09:22 2013 -0500

    BESSEL[IJKY]: Improve accuracy.
    
    Basically replace sin(x*Pi) with sin(fmod(x,2)*Pi) because the argument
    reduction in the latter case can be done with no loss of accuracy.

 ChangeLog      |    5 +++++
 NEWS           |    1 +
 src/mathfunc.c |   18 +++++++++---------
 3 files changed, 15 insertions(+), 9 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 9aae262..64ce20d 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2013-11-03  Morten Welinder  <terra gnome org>
+
+       * src/mathfunc.c (bessel_i, etc.): Do argument reduction for
+       sin/cos before scaling by pi.
+
 2013-11-01  Morten Welinder  <terra gnome org>
 
        * src/mathfunc.c (dpois_raw): Handler x=0 as in newer R.
diff --git a/NEWS b/NEWS
index 6c793b7..569f9d9 100644
--- a/NEWS
+++ b/NEWS
@@ -8,6 +8,7 @@ Morten:
        * Improve accuracy of R.QGAMMA and thus R.QCHISQ.
        * Improve accuracy of R.QBETA, R.QF, R.QTUKEY, R.QSNORM, and R.QST.
        * Improve accuracy of COMBIN, PERMUT, POCHHAMMER, FACT, GAMMA.
+       * Improve accuracy of bessel functions with large non-integer alpha.
 
 Xabier Rodríguez Calvar:
        * Fix dialog button order. [#710378]
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 1d91ce1..3cb8cb7 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -3887,7 +3887,7 @@ static gnm_float bessel_i(gnm_float x, gnm_float alpha, gnm_float expo)
         * this may not be quite optimal (CPU and accuracy wise) */
        return(bessel_i(x, -alpha, expo) +
               bessel_k(x, -alpha, expo) * ((ize == 1)? 2. : 2.*gnm_exp(-x))/M_PIgnum
-              * gnm_sin(-M_PIgnum * alpha));
+              * gnm_sin(-M_PIgnum * gnm_fmod (alpha, 2)));
     }
     nb = 1+ (long)gnm_floor(alpha);/* nb-1 <= alpha < nb */
     alpha -= (nb-1);
@@ -4362,9 +4362,9 @@ static gnm_float bessel_j(gnm_float x, gnm_float alpha)
     if (alpha < 0) {
        /* Using Abramowitz & Stegun  9.1.2
         * this may not be quite optimal (CPU and accuracy wise) */
-       return(bessel_j(x, -alpha) * gnm_cos(M_PIgnum * alpha) +
+       return(bessel_j(x, -alpha) * gnm_cos(M_PIgnum * gnm_fmod (alpha, 2)) +
               ((alpha == na) ? 0 :
-              bessel_y(x, -alpha) * gnm_sin(M_PIgnum * alpha)));
+              bessel_y(x, -alpha) * gnm_sin(M_PIgnum * gnm_fmod (alpha, 2))));
     }
     nb = 1 + (long)na; /* nb-1 <= alpha < nb */
     alpha -= (gnm_float)(nb-1);
@@ -4412,9 +4412,9 @@ static gnm_float bessel_j_ex(gnm_float x, gnm_float alpha, gnm_float *bj)
     if (alpha < 0) {
        /* Using Abramowitz & Stegun  9.1.2
         * this may not be quite optimal (CPU and accuracy wise) */
-       return(bessel_j_ex(x, -alpha, bj) * gnm_cos(M_PIgnum * alpha) +
+       return(bessel_j_ex(x, -alpha, bj) * gnm_cos(M_PIgnum * gnm_fmod (alpha, 2)) +
               ((alpha == na) ? 0 :
-               bessel_y_ex(x, -alpha, bj) * gnm_sin(M_PIgnum * alpha)));
+               bessel_y_ex(x, -alpha, bj) * gnm_sin(M_PIgnum * gnm_fmod (alpha, 2))));
     }
     nb = 1 + (long)na; /* nb-1 <= alpha < nb */
     alpha -= (gnm_float)(nb-1);
@@ -5462,9 +5462,9 @@ static gnm_float bessel_y(gnm_float x, gnm_float alpha)
     if (alpha < 0) {
        /* Using Abramowitz & Stegun  9.1.2
         * this may not be quite optimal (CPU and accuracy wise) */
-       return(bessel_y(x, -alpha) * gnm_cos(M_PIgnum * alpha) -
+       return(bessel_y(x, -alpha) * gnm_cos(M_PIgnum * gnm_fmod (alpha, 2)) -
               ((alpha == na) ? 0 :
-               bessel_j(x, -alpha) * gnm_sin(M_PIgnum * alpha)));
+               bessel_j(x, -alpha) * gnm_sin(M_PIgnum * gnm_fmod (alpha, 2))));
     }
     nb = 1+ (long)na;/* nb-1 <= alpha < nb */
     alpha -= (gnm_float)(nb-1);
@@ -5514,9 +5514,9 @@ static gnm_float bessel_y_ex(gnm_float x, gnm_float alpha, gnm_float *by)
     if (alpha < 0) {
        /* Using Abramowitz & Stegun  9.1.2
         * this may not be quite optimal (CPU and accuracy wise) */
-       return(bessel_y_ex(x, -alpha, by) * gnm_cos(M_PIgnum * alpha) -
+       return(bessel_y_ex(x, -alpha, by) * gnm_cos(M_PIgnum * gnm_fmod (alpha, 2)) -
               ((alpha == na) ? 0 :
-               bessel_j_ex(x, -alpha, by) * gnm_sin(M_PIgnum * alpha)));
+               bessel_j_ex(x, -alpha, by) * gnm_sin(M_PIgnum * gnm_fmod (alpha, 2))));
     }
     nb = 1+ (long)na;/* nb-1 <= alpha < nb */
     alpha -= (gnm_float)(nb-1);