[gnumeric] POCHHAMMER: Improve precision a bit.

[Morten Welinder <mortenw src gnome org>]

commit b9cc2ed4f1187f293e695ad51f0f1039ac38f970
Author: Morten Welinder <terra gnome org>
Date:   Sat Apr 13 14:17:18 2013 -0400

    POCHHAMMER: Improve precision a bit.
    
    This avoids the cancellation involved in (simplified)
    
        x log(x+n) - x log(x)
    
    Note, that this cancellation would grow with x.

 ChangeLog      |  2 ++
 src/mathfunc.c | 23 ++++++++++++-----------
 2 files changed, 14 insertions(+), 11 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 34c3981..80dbd02 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,7 @@
 2013-04-13  Morten Welinder  <terra gnome org>
 
+       * src/mathfunc.c (pochhammer): Improve precision a bit.
+
        * configure.ac (GETTEXT_PACKAGE): Add version number so that
        doesn't prevent multiple versions from co-existing.
 
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 7f52322..a14e1ab 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -7226,19 +7226,20 @@ pochhammer (gnm_float x, gnm_float n, gboolean give_log)
        }
 
        if (give_log) {
-               gnm_float lr;
+               gnm_float lf = 0, lr;
 
-               if (x < 10 || x + n < 10) {
-                       int s;
-                       lr = gnm_lgamma_r (x + n, &s) - gnm_lgamma_r (x, &s);
-               } else {
-                       lr = ((x + n - 0.5) * gnm_log (x + n) -
-                             (x     - 0.5) * gnm_log (x) -
-                             n +
-                             lgammacor (x + n) -
-                             lgammacor (x));
+               while (x < 10 || x + n < 10) {
+                       /* P(x,n) = (x/(x+n)) * P(x+1,n) */
+                       lf -= gnm_log1p (n / x);
+                       x++;
                }
-               return lr;
+
+               lr = ((x - 0.5) * gnm_log1p (n / x) +
+                     n * gnm_log (x + n) -
+                     n +
+                     (lgammacor (x + n) - lgammacor (x)));
+
+               return lr + lf;
        } else {
                gnm_float lr = pochhammer (x, n, TRUE);
                return gnm_exp (lr);