[gnumeric] fn-r: implement a bit more of the skewed t distribution.

[Morten Welinder <mortenw src gnome org>]

commit b30b7e92476ab5bbe893b7d904c4496901cf198d
Author: Morten Welinder <terra gnome org>
Date:   Wed Feb 27 16:48:46 2013 -0500

    fn-r: implement a bit more of the skewed t distribution.

 plugins/fn-r/ChangeLog |    5 +++
 plugins/fn-r/extra.c   |   83 +++++++++++++++++++++++++++++++++++++----------
 2 files changed, 70 insertions(+), 18 deletions(-)
---
diff --git a/plugins/fn-r/ChangeLog b/plugins/fn-r/ChangeLog
index fdac486..73a9eb6 100644
--- a/plugins/fn-r/ChangeLog
+++ b/plugins/fn-r/ChangeLog
@@ -1,3 +1,8 @@
+2013-02-27  Morten Welinder  <terra gnome org>
+
+       * extra.c (qsnorm): Implement.
+       (qst): Tentatively implement based on pst.
+
 2012-12-18  Morten Welinder <terra gnome org>
 
        * Release 1.12.0
diff --git a/plugins/fn-r/extra.c b/plugins/fn-r/extra.c
index ffc28ce..758a945 100644
--- a/plugins/fn-r/extra.c
+++ b/plugins/fn-r/extra.c
@@ -137,7 +137,7 @@ dsnorm (gnm_float x, gnm_float shape, gnm_float location, gnm_float scale, gbool
        if (shape == 0.)
                return dnorm (x, location, scale, give_log);
        else if (give_log)
-               return gnm_log (2.) + dnorm (x, location, scale, TRUE) + pnorm (shape * x, shape * location, 
scale, TRUE, TRUE);
+               return M_LN2gnum + dnorm (x, location, scale, TRUE) + pnorm (shape * x, shape * location, 
scale, TRUE, TRUE);
        else
                return 2 * dnorm (x, location, scale, FALSE) * pnorm (shape * x, location/shape, scale, TRUE, 
FALSE);
 }
@@ -161,16 +161,36 @@ psnorm (gnm_float x, gnm_float shape, gnm_float location, gnm_float scale, gbool
                return result;
 }
 
+static gnm_float
+dsnorm1 (gnm_float x, const gnm_float params[], gboolean log_p)
+{
+       return dsnorm (x, params[0], params[1], params[2], log_p);
+}
+
+static gnm_float
+psnorm1 (gnm_float x, const gnm_float params[],
+        gboolean lower_tail, gboolean log_p)
+{
+       return psnorm (x, params[0], params[1], params[2], lower_tail, log_p);
+}
 
 gnm_float
-qsnorm (gnm_float p, gnm_float shape, gnm_float location, gnm_float scale, gboolean lower_tail, gboolean 
log_p)
+qsnorm (gnm_float p, gnm_float shape, gnm_float location, gnm_float scale,
+       gboolean lower_tail, gboolean log_p)
 {
+       gnm_float x0;
+       gnm_float params[3];
+
        if (shape == 0.)
                return qnorm (p, location, scale, lower_tail, log_p);
-       else if (log_p)
-               return 0.;
-       else
-               return 0;
+
+       x0 = 0.0;
+       params[0] = shape;
+       params[1] = location;
+       params[2] = scale;
+       return pfuncinverter (p, params, lower_tail, log_p,
+                             gnm_ninf, gnm_pinf, x0,
+                             psnorm1, dsnorm1);
 }
 
 /* ------------------------------------------------------------------------- */
@@ -186,7 +206,7 @@ dst (gnm_float x, gnm_float n, gnm_float shape, gboolean give_log)
                gnm_float pdf = dt (x, n, give_log);
                gnm_float cdf = pt (shape * x * gnm_sqrt ((n + 1)/(x * x + n)),
                                    n + 1, TRUE, give_log);
-               return ((give_log) ? (gnm_log (2.) + pdf + cdf) : (2. * pdf * cdf));
+               return give_log ? (M_LN2gnum + pdf + cdf) : (2. * pdf * cdf);
        }
 }
 
@@ -196,25 +216,52 @@ pst (gnm_float x, gnm_float n, gnm_float shape, gboolean lower_tail, gboolean lo
 {
        if (shape == 0.)
                return pt (x, n, lower_tail, log_p);
-       else if (log_p)
-               return 0.;
-       else
-               return 0.;
+
+       /* Generic fallback.  */
+       if (!lower_tail)
+               return log_p
+                       ? swap_log_tail (pst (x, n, shape, TRUE, TRUE))
+                       : 1 - pst (x, n, shape, TRUE, FALSE);
+       if (log_p)
+               gnm_log (pst (x, n, shape, TRUE, FALSE));
+
+       /* FIXME: Numerical integration of dst from -inf to x.  */
+
+       /* Give up.  */
+       return gnm_nan;
 }
 
 
+static gnm_float
+dst1 (gnm_float x, const gnm_float params[], gboolean log_p)
+{
+       return dst (x, params[0], params[1], log_p);
+}
+
+static gnm_float
+pst1 (gnm_float x, const gnm_float params[],
+      gboolean lower_tail, gboolean log_p)
+{
+       return pst (x, params[0], params[1], lower_tail, log_p);
+}
+
 gnm_float
-qst (gnm_float p, gnm_float n, gnm_float shape, gboolean lower_tail, gboolean log_p)
+qst (gnm_float p, gnm_float n, gnm_float shape,
+     gboolean lower_tail, gboolean log_p)
 {
+       gnm_float x0;
+       gnm_float params[2];
+
        if (shape == 0.)
                return qt (p, n, lower_tail, log_p);
-       else if (log_p)
-               return 0.;
-       else
-               return 0.;
-}
-
 
+       x0 = 0.0;
+       params[0] = n;
+       params[1] = shape;
+       return pfuncinverter (p, params, lower_tail, log_p,
+                             gnm_ninf, gnm_pinf, x0,
+                             pst1, dst1);
+}
 
 /* ------------------------------------------------------------------------- */