[gnumeric] dpois: improve accuracy.

[Morten Welinder <mortenw src gnome org>]

commit 18553c96715c0fe5f4f8e71a5f7899a000176103
Author: Morten Welinder <terra gnome org>
Date:   Mon Dec 30 16:49:23 2013 -0500

    dpois: improve accuracy.

 ChangeLog              |    2 +
 samples/amath.gnumeric |  Bin 44226 -> 44337 bytes
 src/mathfunc.c         |   95 ++++++++++++++++++------------------------------
 3 files changed, 37 insertions(+), 60 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 98a200d..5309379 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,6 +1,8 @@
 2013-12-30  Morten Welinder  <terra gnome org>
 
        * src/mathfunc.c (dnorm): Improve accuracy for x>5 (normalized).
+       (bd0): Reimplement.
+       (dpois_raw): Avoid going through logs, if possible.
 
 2013-12-25  Morten Welinder  <terra gnome org>
 
diff --git a/samples/amath.gnumeric b/samples/amath.gnumeric
index 940e18e..9ac54ba 100644
Binary files a/samples/amath.gnumeric and b/samples/amath.gnumeric differ
diff --git a/src/mathfunc.c b/src/mathfunc.c
index aea76f9..ade4892 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -77,6 +77,15 @@ mathfunc_init (void)
        /* Nothing, for the time being.  */
 }
 
+static gnm_float
+bd0(gnm_float x, gnm_float M)
+{
+       return (gnm_abs (M - x) < x / 4)
+               ? -x * log1pmx ((M - x) / x)
+               : x * gnm_log (x / M) + (M - x);
+}
+
+
 /* ------------------------------------------------------------------------- */
 /* --- BEGIN MAGIC R SOURCE MARKER --- */
 
@@ -826,66 +835,6 @@ gnm_float ppois(gnm_float x, gnm_float lambda, gboolean lower_tail, gboolean log
 }
 
 /* ------------------------------------------------------------------------ */
-/* Imported src/nmath/bd0.c from R.  */
-/*
- *  AUTHOR
- *     Catherine Loader, catherine research bell-labs com 
- *     October 23, 2000.
- *
- *  Merge in to R:
- *     Copyright (C) 2000, The R Core Development Team
- *
- *  This program is free software; you can redistribute it and/or modify
- *  it under the terms of the GNU General Public License as published by
- *  the Free Software Foundation; either version 2 of the License, or
- *  (at your option) any later version.
- *
- *  This program is distributed in the hope that it will be useful,
- *  but WITHOUT ANY WARRANTY; without even the implied warranty of
- *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- *  GNU General Public License for more details.
- *
- *  You should have received a copy of the GNU General Public License
- *  along with this program; if not, write to the Free Software
- *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
- *  USA.
- *
- *
- *  DESCRIPTION
- *     Evaluates the "deviance part"
- *     bd0(x,M) :=  M * D0(x/M) = M*[ x/M * log(x/M) + 1 - (x/M) ] =
- *               =  x * log(x/M) + M - x
- *     where M = E[X] = n*p (or = lambda), for   x, M > 0
- *
- *     in a manner that should be stable (with small relative error)
- *     for all x and np. In particular for x/np close to 1, direct
- *     evaluation fails, and evaluation is based on the Taylor series
- *     of log((1+v)/(1-v)) with v = (x-np)/(x+np).
- */
-
-static gnm_float bd0(gnm_float x, gnm_float np)
-{
-    gnm_float ej, s, s1, v;
-    int j;
-
-    if (gnm_abs(x-np) < 0.1*(x+np)) {
-       v = (x-np)/(x+np);
-       s = (x-np)*v;/* s using v -- change by MM */
-       ej = 2*x*v;
-       v = v*v;
-       for (j=1; ; j++) { /* Taylor series */
-           ej *= v;
-           s1 = s+ej/((j<<1)+1);
-           if (s1==s) /* last term was effectively 0 */
-               return(s1);
-           s = s1;
-       }
-    }
-    /* else:  | x - np |  is not too small */
-    return(x*gnm_log(x/np)+np-x);
-}
-
-/* ------------------------------------------------------------------------ */
 /* Imported src/nmath/dpois.c from R.  */
 /*
  *  AUTHOR
@@ -925,6 +874,9 @@ static gnm_float bd0(gnm_float x, gnm_float np)
 
 static gnm_float dpois_raw(gnm_float x, gnm_float lambda, gboolean give_log)
 {
+    GnmQuad mfx;
+    int efx;
+
     /*       x >= 0 ; integer for dpois(), but not e.g. for pgamma()!
         lambda >= 0
     */
@@ -933,6 +885,29 @@ static gnm_float dpois_raw(gnm_float x, gnm_float lambda, gboolean give_log)
     if (x < 0) return( R_D__0 );
     if (x <= lambda * GNM_MIN) return(R_D_exp(-lambda) );
     if (lambda < x * GNM_MIN) return(R_D_exp(-lambda + x*gnm_log(lambda) -lgamma1p (x)));
+
+    if (!qfactf (x, &mfx, &efx)) {
+           void *state = gnm_quad_start ();
+           gnm_float r, elx, eel;
+           GnmQuad qr, qx, ql, mlx, mel;
+
+           gnm_quad_init (&ql, lambda);
+           gnm_quad_init (&qx, x);
+           gnm_quad_pow (&mlx, &elx, &ql, &qx);
+           gnm_quad_exp (&mel, &eel, &ql);
+           gnm_quad_mul (&qr, &mfx, &mel);
+           gnm_quad_div (&qr, &mlx, &qr);
+           r = gnm_quad_value (&qr);
+           gnm_quad_end (state);
+
+           if (gnm_finite (r)) {
+                   gnm_float e = elx - eel - efx;
+                   return give_log
+                           ? gnm_log (r) + M_LN2gnum * e
+                           : gnm_ldexp (r, CLAMP (e, G_MININT, G_MAXINT));
+           }
+    }
+
     return(R_D_fexp( M_2PIgnum*x, -stirlerr(x)-bd0(x,lambda) ));
 }