[gnumeric] Math: detect underflow in direct compilation of dpois.

[Morten Welinder <mortenw src gnome org>]

commit 51cf04ac55aee9c9ad351625e83d49bcaa6ebad3
Author: Morten Welinder <terra gnome org>
Date:   Wed Jan 8 10:56:53 2014 -0500

    Math: detect underflow in direct compilation of dpois.

 ChangeLog      |    2 ++
 src/mathfunc.c |   24 +++++++++++++++++-------
 src/numbers.h  |    2 ++
 3 files changed, 21 insertions(+), 7 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 63e886a..3b60ec3 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,7 @@
 2014-01-08  Morten Welinder  <terra gnome org>
 
+       * src/mathfunc.c (dpois_raw): Detect underflow in direct formula.
+
        * src/gui-file.h: Namespace improvements.  Also make gui_file_read
        return the WorkbookView instead of a boolean.
 
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 4e97fba..d820242 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -743,14 +743,24 @@ static gnm_float dpois_raw(gnm_float x, gnm_float lambda, gboolean give_log)
            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)) {
+           if (gnm_finite (r) && r > 0) {
                    gnm_float e = elx - eel - efx;
-                   return give_log
-                           ? gnm_log (r) + M_LN2gnum * e
-                           : gnm_ldexp (r, CLAMP (e, G_MININT, G_MAXINT));
-           }
+                   if (give_log) {
+                           GnmQuad qt;
+                           gnm_quad_init (&qt, e);
+                           gnm_quad_mul (&qt, &qt, &gnm_quad_ln2);
+                           gnm_quad_log (&qr, &qr);
+                           gnm_quad_add (&qr, &qr, &qt);
+                           r = gnm_quad_value (&qr);
+                   } else {
+                           r = gnm_ldexp (r, CLAMP (e, G_MININT, G_MAXINT));
+                   }
+           } else
+                   r = gnm_ninf;
+           gnm_quad_end (state);
+
+           if (gnm_finite (r))
+                   return r;
     }
 
     return(R_D_fexp( M_2PIgnum*x, -stirlerr(x)-bd0(x,lambda) ));
diff --git a/src/numbers.h b/src/numbers.h
index cc9fb11..9c374e1 100644
--- a/src/numbers.h
+++ b/src/numbers.h
@@ -115,6 +115,7 @@ typedef long double gnm_float;
 #define gnm_quad_exp go_quad_expl
 #define gnm_quad_expm1 go_quad_expm1l
 #define gnm_quad_init go_quad_initl
+#define gnm_quad_ln2 go_quad_ln2l
 #define gnm_quad_log go_quad_logl
 #define gnm_quad_mul go_quad_mull
 #define gnm_quad_mul12 go_quad_mul12l
@@ -225,6 +226,7 @@ typedef double gnm_float;
 #define gnm_quad_exp go_quad_exp
 #define gnm_quad_expm1 go_quad_expm1
 #define gnm_quad_init go_quad_init
+#define gnm_quad_ln2 go_quad_ln2
 #define gnm_quad_log go_quad_log
 #define gnm_quad_mul go_quad_mul
 #define gnm_quad_mul12 go_quad_mul12