[gnumeric] BETA: Improve accuracy.

[Morten Welinder <mortenw src gnome org>]

commit e6c4e52f1c21c64c177e65633d9812bb8330cce6
Author: Morten Welinder <terra gnome org>
Date:   Fri Jan 29 21:00:37 2016 -0500

    BETA: Improve accuracy.
    
    When one argument is less than one, such as in BETA(1000,1e-10), it's
    better to use pochhammer because forming the sum is problematic.

 ChangeLog      |    3 +++
 NEWS           |    1 +
 src/sf-gamma.c |   19 +++++++++++++++++++
 3 files changed, 23 insertions(+), 0 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 5c85c66..14a1d55 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,8 @@
 2016-01-29  Morten Welinder  <terra gnome org>
 
+       * src/sf-gamma.c (qbetaf): Improve accuracy in the case where one
+       argument is less than one.
+
        * src/xml-sax-read.c (xml_sax_filter_condition): Leak fix and warn
        about broken sheet filter.
 
diff --git a/NEWS b/NEWS
index dc8e17a..af291cc 100644
--- a/NEWS
+++ b/NEWS
@@ -23,6 +23,7 @@ Morten:
        * Fix canvas problem leaving grab in place.  [#760639]
        * Work around gtk+ bug causing growing windows.  [#761142]
        * Improve BESSELJ and BESSELY.
+       * Improve BETA accuracy.
 
 --------------------------------------------------------------------------
 Gnumeric 1.12.26
diff --git a/src/sf-gamma.c b/src/sf-gamma.c
index cbebbcd..ae46197 100644
--- a/src/sf-gamma.c
+++ b/src/sf-gamma.c
@@ -8,6 +8,7 @@
 #define ML_ERROR(cause) do { } while(0)
 
 static int qgammaf (gnm_float x, GnmQuad *mant, int *exp2);
+static void pochhammer_small_n (gnm_float x, gnm_float n, GnmQuad *res);
 
 
 /* Compute  gnm_log(gamma(a+1))  accurately also for small a (0 < a < 0.5). */
@@ -515,6 +516,24 @@ qbetaf (gnm_float a, gnm_float b, GnmQuad *mant, int *exp2)
                return 0;
        }
 
+       if (b > a) {
+               gnm_float s = a;
+               a = b;
+               b = s;
+       }
+
+       if (a > 20 && gnm_abs (b) < 1) {
+               void *state;
+               if (qgammaf (b, &mb, &eb))
+                       return 1;
+               state = gnm_quad_start ();
+               pochhammer_small_n (a, b, &ma);
+               gnm_quad_div (mant, &mb, &ma);
+               gnm_quad_end (state);
+               *exp2 = eb;
+               return 0;
+       }
+
        if (!qgammaf (a, &ma, &ea) &&
            !qgammaf (b, &mb, &eb) &&
            !qgammaf (ab, &mab, &eab)) {