[gnumeric] pnorm2: drop log_p argument.

[Morten Welinder <mortenw src gnome org>]

commit bd539943fe51bc012eb30590db5bbb522c4a2737
Author: Morten Welinder <terra gnome org>
Date:   Mon May 27 18:34:40 2013 -0400

    pnorm2: drop log_p argument.

 ChangeLog                |    5 ++++
 plugins/fn-r/functions.c |    6 +---
 src/mathfunc.c           |   53 ++++++++++++++++++---------------------------
 src/mathfunc.h           |    2 +-
 4 files changed, 29 insertions(+), 37 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index b4a126c..30c94c2 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2013-05-27  Morten Welinder  <terra gnome org>
+
+       * src/mathfunc.c (pnorm2): Drop log_p argument.  We don't need it
+       and we aren't very good at it.
+
 2013-05-27  Andreas J. Guelzow <aguelzow pyrshep ca>
 
        * component/Gnumeric-embed.xml.in: ctrl-a will be handled in gnm-pane.c
diff --git a/plugins/fn-r/functions.c b/plugins/fn-r/functions.c
index d5f20c8..c509178 100644
--- a/plugins/fn-r/functions.c
+++ b/plugins/fn-r/functions.c
@@ -71,7 +71,6 @@ static GnmFuncHelp const help_r_pnorm2[] = {
        { GNM_FUNC_HELP_NAME, F_("R.PNORM2:cumulative distribution function of the normal distribution") },
        { GNM_FUNC_HELP_ARG, F_("x1:first observation") },
        { GNM_FUNC_HELP_ARG, F_("x2:first observation") },
-       { GNM_FUNC_HELP_ARG, F_("log_p:if true, log of the probability is used") },
        { GNM_FUNC_HELP_DESCRIPTION, F_("This function returns the cumulative distribution function of the 
normal distribution.") },
        { GNM_FUNC_HELP_END }
 };
@@ -81,9 +80,8 @@ gnumeric_r_pnorm2 (GnmFuncEvalInfo *ei, GnmValue const * const *args)
 {
        gnm_float x1 = value_get_as_float (args[0]);
        gnm_float x2 = value_get_as_float (args[1]);
-       gboolean log_p = args[2] ? value_get_as_checked_bool (args[2]) : FALSE;
 
-       return value_new_float (pnorm2 (x1, x2, log_p));
+       return value_new_float (pnorm2 (x1, x2));
 }
 
 /* ------------------------------------------------------------------------- */
@@ -1451,7 +1449,7 @@ GnmFuncDescriptor const stat_functions[] = {
        },
        {
                "r.pnorm2",
-               "ff|b",
+               "ff",
                help_r_pnorm2,
                gnumeric_r_pnorm2, NULL, NULL, NULL,
                GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE,
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 78bd027..aa4bcad 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -5243,7 +5243,7 @@ static gnm_float ptukey_wprob(gnm_float w, gnm_float rr, gnm_float cc)
            }
            v = C + binc * 0.5 * xx;
 
-           rinsum = pnorm2 (v - w, v, FALSE);
+           rinsum = pnorm2 (v - w, v);
            elsum += gnm_pow (rinsum, cc - 1) *
                    aa *
                    gnm_exp(-0.5 * v * v);
@@ -5589,61 +5589,50 @@ swap_log_tail (gnm_float lp)
 
 
 gnm_float
-pnorm2 (gnm_float x1, gnm_float x2, gboolean log_p)
+pnorm2 (gnm_float x1, gnm_float x2)
 {
        if (gnm_isnan(x1) || gnm_isnan(x2))
                return gnm_nan;
 
        if (x1 > x2)
-               return log_p ? gnm_nan : 0 - pnorm2 (x2, x1, log_p);
+               return 0 - pnorm2 (x2, x1);
 
+       /* A bunch of special cases:  */
        if (x1 == x2)
-               return R_D__0;
-
+               return 0.0;
        if (x1 == gnm_ninf)
-               return pnorm (x2, 0.0, 1.0, TRUE, log_p);
-
+               return pnorm (x2, 0.0, 1.0, TRUE, FALSE);
        if (x2 == gnm_pinf)
-               return pnorm (x1, 0.0, 1.0, FALSE, log_p);
-
-       if (x1 == 0 && !log_p) {
+               return pnorm (x1, 0.0, 1.0, FALSE, FALSE);
+       if (x1 == 0)
                return gnm_erf (x2 / M_SQRT2gnum) / 2;
-
-               /*
-                * for the log case we can use this log(erf(x)) expansion:
-                * (log(2 x)-(log(pi))/2)-x^2/3+(2 x^4)/45-(8 x^6)/2835-(4 x^8)/14175+(32 x^10)/467775+(736 
x^12)/1915538625-(2944 x^14)/1915538625+(5024 x^16)/44405668125+(49690112 x^18)/1754068296605625+O(x^20)
-                */
-       }
-       if (x2 == 0 && !log_p) {
+       if (x2 == 0)
                return gnm_erf (x1 / -M_SQRT2gnum) / 2;
-       }
 
        if (x1 <= 0 && x2 >= 0) {
                /* The interval spans 0.  */
-               gnm_float p1 = pnorm2 (0, MIN (-x1, x2), log_p);
-               gnm_float p2 = pnorm2 (MIN (-x1, x2), MAX (-x1, x2), log_p);
-               return log_p
-                       ? logspace_add (p1 + M_LN2gnum, p2)
-                       : 2 * p1 + p2;
+               gnm_float p1 = pnorm2 (0, MIN (-x1, x2));
+               gnm_float p2 = pnorm2 (MIN (-x1, x2), MAX (-x1, x2));
+               return 2 * p1 + p2;
        } else if (x1 < 0) {
                /* Both < 0 -- use symmetry */
-               return pnorm2 (-x2, -x1, log_p);
+               return pnorm2 (-x2, -x1);
        } else {
                /* Both >= 0 */
-               gnm_float p1C = pnorm (x1, 0.0, 1.0, FALSE, log_p);
-               gnm_float p2C = pnorm (x2, 0.0, 1.0, FALSE, log_p);
-               gnm_float raw = log_p ? logspace_sub (p1C, p2C) : p1C - p2C;
+               gnm_float p1C = pnorm (x1, 0.0, 1.0, FALSE, FALSE);
+               gnm_float p2C = pnorm (x2, 0.0, 1.0, FALSE, FALSE);
+               gnm_float raw = p1C - p2C;
                gnm_float dx, d1, d2, ub, lb;
 
                if (gnm_abs (p1C - p2C) * 32 > gnm_abs (p1C + p2C))
                        return raw;
 
                /* dnorm is strictly decreasing in this area.  */
-               dx = log_p ? gnm_log (x2 - x1) : x2 - x1;
-               d1 = dnorm (x1, 0.0, 1.0, log_p);
-               d2 = dnorm (x2, 0.0, 1.0, log_p);
-               ub = log_p ? dx + d1 : dx * d1;  /* upper bound */
-               lb = log_p ? dx + d2 : dx * d2;  /* lower bound */
+               dx = x2 - x1;
+               d1 = dnorm (x1, 0.0, 1.0, FALSE);
+               d2 = dnorm (x2, 0.0, 1.0, FALSE);
+               ub = dx * d1;  /* upper bound */
+               lb = dx * d2;  /* lower bound */
 
                raw = MAX (raw, lb);
                raw = MIN (raw, ub);
diff --git a/src/mathfunc.h b/src/mathfunc.h
index d920009..2c15b61 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -56,7 +56,7 @@ gnm_float bessel_k (gnm_float x, gnm_float alpha, gnm_float expo);
 /* The normal distribution.  */
 gnm_float dnorm (gnm_float x, gnm_float mu, gnm_float sigma, gboolean give_log);
 gnm_float pnorm (gnm_float x, gnm_float mu, gnm_float sigma, gboolean lower_tail, gboolean log_p);
-gnm_float pnorm2 (gnm_float x1, gnm_float x2, gboolean log_p);
+gnm_float pnorm2 (gnm_float x1, gnm_float x2);
 gnm_float qnorm (gnm_float p, gnm_float mu, gnm_float sigma, gboolean lower_tail, gboolean log_p);
 
 /* The log-normal distribution.  */