[gnumeric] BESSELJ, BESSELY: Handle fractional order.

[Morten Welinder <mortenw src gnome org>]

commit 9b8df6bf5793ce1071db1651703884fc5b559c6b
Author: Morten Welinder <terra gnome org>
Date:   Wed Aug 28 14:26:01 2013 -0400

    BESSELJ, BESSELY: Handle fractional order.
    
    This is an extension compared to XL.

 NEWS                       |    3 ++
 plugins/fn-eng/functions.c |   51 ++++++++++---------------------------------
 src/mathfunc.c             |   52 ++++++++++++++++++++++++++++++++++++++++---
 src/mathfunc.h             |    8 +++---
 4 files changed, 67 insertions(+), 47 deletions(-)
---
diff --git a/NEWS b/NEWS
index 2458fd8..17e12be 100644
--- a/NEWS
+++ b/NEWS
@@ -1,5 +1,8 @@
 Gnumeric 1.12.7
 
+Morten:
+       * Extend BESSELJ and BESSELY to handle fractional order.  [#706720]
+
 --------------------------------------------------------------------------
 Gnumeric 1.12.6
 
diff --git a/plugins/fn-eng/functions.c b/plugins/fn-eng/functions.c
index 3763a94..ab13903 100644
--- a/plugins/fn-eng/functions.c
+++ b/plugins/fn-eng/functions.c
@@ -555,25 +555,9 @@ static GnmFuncHelp const help_besseli[] = {
 static GnmValue *
 gnumeric_besseli (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
 {
-       gnm_float x = value_get_as_float (argv[0]);     /* value to evaluate I_n at. */
-       gnm_float order = value_get_as_float (argv[1]); /* the order */
-       gnm_float r;
-
-       if (order < 0)
-               return value_new_error_NUM (ei->pos);
-
-       /* This, or something like it, ought to be moved into a proper bessel_i.  */
-       if (x < 0) {
-               if (order != gnm_floor (order))
-                       return value_new_error_NUM (ei->pos);
-               else if (gnm_fmod (order, 2) == 0)
-                       r = bessel_i (-x, order, 1);  /* Even for even order */
-               else
-                       r = -bessel_i (-x, order, 1);  /* Odd for odd order */
-       } else
-               r = bessel_i (x, order, 1);
-
-       return value_new_float (r);
+       gnm_float x = value_get_as_float (argv[0]);
+       gnm_float order = value_get_as_float (argv[1]);
+       return value_new_float (gnm_bessel_i (x, order));
 }
 
 /***************************************************************************/
@@ -594,10 +578,9 @@ static GnmFuncHelp const help_besselk[] = {
 static GnmValue *
 gnumeric_besselk (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
 {
-       gnm_float x = value_get_as_float (argv[0]);     /* value to evaluate K_n at. */
+       gnm_float x = value_get_as_float (argv[0]);
        gnm_float order = value_get_as_float (argv[1]); /* the order */
-
-       return value_new_float (bessel_k (x, order, 1.0));
+       return value_new_float (gnm_bessel_k (x, order));
 }
 
 /***************************************************************************/
@@ -607,9 +590,8 @@ static GnmFuncHelp const help_besselj[] = {
         { GNM_FUNC_HELP_ARG, F_("X:number") },
         { GNM_FUNC_HELP_ARG, F_("\xce\xb1:order (any non-negative integer)") },
        { GNM_FUNC_HELP_NOTE, F_("If @{x} or @{\xce\xb1} are not numeric, #VALUE! is returned. "
-                                "If @{\xce\xb1} < 0, #NUM! is returned. "
-                                "If @{\xce\xb1} is not an integer, it is truncated.") },
-       { GNM_FUNC_HELP_EXCEL, F_("This function is Excel compatible.") },
+                                "If @{\xce\xb1} < 0, #NUM! is returned.") },
+       { GNM_FUNC_HELP_EXCEL, F_("This function is Excel compatible if only integer orders @{\xce\xb1} are 
used.") },
         { GNM_FUNC_HELP_EXAMPLES, "=BESSELJ(0.89,3)" },
         { GNM_FUNC_HELP_SEEALSO, "BESSELI,BESSELK,BESSELY" },
        { GNM_FUNC_HELP_EXTREF, F_("wiki:en:Bessel_function") },
@@ -621,12 +603,7 @@ gnumeric_besselj (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
 {
        gnm_float x = value_get_as_float (argv[0]);
        gnm_float y = value_get_as_float (argv[1]);
-
-       if (y < 0 || y > INT_MAX)
-               return value_new_error_NUM (ei->pos);
-
-       /* FIXME: Why not gnm_jn?  */
-       return value_new_float (jn ((int)y, x));
+       return value_new_float (gnm_bessel_j (x, y));
 }
 
 /***************************************************************************/
@@ -636,11 +613,11 @@ static GnmFuncHelp const help_bessely[] = {
         { GNM_FUNC_HELP_ARG, F_("X:number") },
         { GNM_FUNC_HELP_ARG, F_("\xce\xb1:order (any non-negative integer)") },
        { GNM_FUNC_HELP_NOTE, F_("If @{x} or @{\xce\xb1} are not numeric, #VALUE! is returned. "
-                                "If @{\xce\xb1} < 0, #NUM! is returned. "
-                                "If @{\xce\xb1} is not an integer, it is truncated.") },
-       { GNM_FUNC_HELP_EXCEL, F_("This function is Excel compatible.") },
+                                "If @{\xce\xb1} < 0, #NUM! is returned.") },
+       { GNM_FUNC_HELP_EXCEL, F_("This function is Excel compatible if only integer orders @{\xce\xb1} are 
used.") },
         { GNM_FUNC_HELP_EXAMPLES, "=BESSELY(4,2)" },
         { GNM_FUNC_HELP_SEEALSO, "BESSELI,BESSELJ,BESSELK" },
+       { GNM_FUNC_HELP_EXTREF, F_("wiki:en:Bessel_function") },
         { GNM_FUNC_HELP_END}
 };
 
@@ -649,11 +626,7 @@ gnumeric_bessely (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
 {
        gnm_float x = value_get_as_float (argv[0]);
        gnm_float y = value_get_as_float (argv[1]);
-
-       if (y < 0 || y > INT_MAX)
-               return value_new_error_NUM (ei->pos);
-
-       return value_new_float (gnm_yn ((int)y, x));
+       return value_new_float (gnm_bessel_y (x, y));
 }
 
 /***************************************************************************/
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 001fbc8..68c8d1a 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -78,6 +78,8 @@ static inline int imax2 (int x, int y) { return MAX (x, y); }
 #define ML_ERR_return_NAN { return gnm_nan; }
 static void pnorm_both (gnm_float x, gnm_float *cum, gnm_float *ccum, int i_tail, gboolean log_p);
 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);
 
 /* MW ---------------------------------------------------------------------- */
 
@@ -3863,7 +3865,7 @@ And routine specific :
 static void I_bessel(gnm_float *x, gnm_float *alpha, long *nb,
                     long *ize, gnm_float *bi, long *ncalc);
 
-gnm_float bessel_i(gnm_float x, gnm_float alpha, gnm_float expo)
+static gnm_float bessel_i(gnm_float x, gnm_float alpha, gnm_float expo)
 {
     long nb, ncalc, ize;
     gnm_float *bi;
@@ -4340,7 +4342,7 @@ L230:
 static void J_bessel(gnm_float *x, gnm_float *alpha, long *nb,
                     gnm_float *b, long *ncalc);
 
-gnm_float bessel_j(gnm_float x, gnm_float alpha)
+static gnm_float bessel_j(gnm_float x, gnm_float alpha)
 {
     long nb, ncalc;
     gnm_float na, *bj;
@@ -4913,7 +4915,7 @@ L250:
 static void K_bessel(gnm_float *x, gnm_float *alpha, long *nb,
                     long *ize, gnm_float *bk, long *ncalc);
 
-gnm_float bessel_k(gnm_float x, gnm_float alpha, gnm_float expo)
+static gnm_float bessel_k(gnm_float x, gnm_float alpha, gnm_float expo)
 {
     long nb, ncalc, ize;
     gnm_float *bk;
@@ -5440,7 +5442,7 @@ L420:
 static void Y_bessel(gnm_float *x, gnm_float *alpha, long *nb,
                     gnm_float *by, long *ncalc);
 
-gnm_float bessel_y(gnm_float x, gnm_float alpha)
+static gnm_float bessel_y(gnm_float x, gnm_float alpha)
 {
     long nb, ncalc;
     gnm_float na, *by;
@@ -8870,3 +8872,45 @@ gnm_gamma (gnm_float x)
 }
 
 /* ------------------------------------------------------------------------- */
+
+gnm_float
+gnm_bessel_i (gnm_float x, gnm_float alpha)
+{
+       if (x < 0) {
+               if (alpha != gnm_floor (alpha))
+                       return gnm_nan;
+
+               return gnm_fmod (alpha, 2) == 0
+                       ? bessel_i (-x, alpha, 1)  /* Even for even alpha */
+                       : 0 - bessel_i (-x, alpha, 1);  /* Odd for odd alpha */
+       } else
+               return bessel_i (x, alpha, 1);
+}
+
+gnm_float
+gnm_bessel_j (gnm_float x, gnm_float alpha)
+{
+       if (x < 0) {
+               if (alpha != gnm_floor (alpha))
+                       return gnm_nan;
+
+               return gnm_fmod (alpha, 2) == 0
+                       ? bessel_j (-x, alpha)  /* Even for even alpha */
+                       : 0 - bessel_j (-x, alpha);  /* Odd for odd alpha */
+       } else
+               return bessel_j (x, alpha);
+}
+
+gnm_float
+gnm_bessel_k (gnm_float x, gnm_float alpha)
+{
+       return bessel_k (x, alpha, 1);
+}
+
+gnm_float
+gnm_bessel_y (gnm_float x, gnm_float alpha)
+{
+       return bessel_y (x, alpha);
+}
+
+/* ------------------------------------------------------------------------- */
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 358deee..290f436 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -47,10 +47,10 @@ gnm_float gnm_acoth (gnm_float x);
 gnm_float beta (gnm_float a, gnm_float b);
 gnm_float lbeta3 (gnm_float a, gnm_float b, int *sign);
 
-gnm_float bessel_i (gnm_float x, gnm_float alpha, gnm_float expo);
-gnm_float bessel_j (gnm_float x, gnm_float alpha);
-gnm_float bessel_k (gnm_float x, gnm_float alpha, gnm_float expo);
-gnm_float bessel_y (gnm_float x, gnm_float alpha);
+gnm_float gnm_bessel_i (gnm_float x, gnm_float alpha);
+gnm_float gnm_bessel_j (gnm_float x, gnm_float alpha);
+gnm_float gnm_bessel_k (gnm_float x, gnm_float alpha);
+gnm_float gnm_bessel_y (gnm_float x, gnm_float alpha);
 
 /* "d": density.  */
 /* "p": distribution function.  */