[gnumeric] R.QCAUCHY: Improve accuracy.

[Morten Welinder <mortenw src gnome org>]

commit a838504d6388e57f9f06207092c96573c52d609b
Author: Morten Welinder <terra gnome org>
Date:   Sat Apr 11 22:54:50 2015 -0400

    R.QCAUCHY: Improve accuracy.
    
    With non-zero location we can suffer cancellation.

 ChangeLog      |    4 ++++
 NEWS           |    1 +
 src/mathfunc.c |   12 +++++-------
 src/numbers.h  |    6 ++++++
 src/sf-dpq.c   |   38 +++++++++++++++++++++++++++++++++-----
 5 files changed, 49 insertions(+), 12 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index a193931..5723c8c 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -5,6 +5,10 @@
 
 2015-04-11  Morten Welinder  <terra gnome org>
 
+       * src/sf-dpq.c (qcauchy): Handle cancellation.
+
+       * src/mathfunc.c (pcauchy): Simplify.
+
        * src/sf-dpq.c (dnorm): Improve accuracy in certain far-tail cases.
        (drayleigh): Import from fn-stat.  Rename.  Improve accuracy.
 
diff --git a/NEWS b/NEWS
index e5fcdd4..11aa0b6 100644
--- a/NEWS
+++ b/NEWS
@@ -33,6 +33,7 @@ Morten:
        * Fix sheet filter problem.
        * Minor R.DNORM improvement.
        * Improve RAYLEIGH accuracy.
+       * Improve R.QCAUCHY accuracy.
 
 --------------------------------------------------------------------------
 Gnumeric 1.12.21
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 54b3b83..50fb273 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -3272,13 +3272,11 @@ gnm_float pcauchy(gnm_float x, gnm_float location, gnm_float scale,
 #endif
     if (!lower_tail)
        x = -x;
-    /* for large x, the standard formula suffers from cancellation.
-     * This is from Morten Welinder thanks to  Ian Smith's  atan(1/x) : */
-    if (gnm_abs(x) > 1) {
-       gnm_float y = atan(1/x) / M_PIgnum;
-       return (x > 0) ? R_D_Clog(y) : R_D_val(-y);
-    } else
-       return R_D_val(0.5 + atan(x) / M_PIgnum);
+
+    if (log_p && x > 0)
+           return R_D_Clog(gnm_atan2pi (1, x));
+    else
+           return R_D_val(gnm_atan2pi (1, -x));
 }
 
 /* ------------------------------------------------------------------------ */
diff --git a/src/numbers.h b/src/numbers.h
index 3ef9565..e38f632 100644
--- a/src/numbers.h
+++ b/src/numbers.h
@@ -41,10 +41,13 @@ typedef long double gnm_float;
 #define gnm_asinh asinhl
 #define gnm_atan atanl
 #define gnm_atan2 atan2l
+#define gnm_atanpi go_atanpil
+#define gnm_atan2pi go_atan2pil
 #define gnm_atanh atanhl
 #define gnm_ceil ceill
 #define gnm_cosh coshl
 #define gnm_cospi go_cospil
+#define gnm_cotpi go_cotpil
 #define gnm_erf erfl
 #define gnm_erfc erfcl
 #define gnm_exp expl
@@ -153,10 +156,13 @@ typedef double gnm_float;
 #define gnm_asinh asinh
 #define gnm_atan atan
 #define gnm_atan2 atan2
+#define gnm_atanpi go_atanpi
+#define gnm_atan2pi go_atan2pi
 #define gnm_atanh atanh
 #define gnm_ceil ceil
 #define gnm_cosh cosh
 #define gnm_cospi go_cospi
+#define gnm_cotpi go_cotpi
 #define gnm_erf erf
 #define gnm_erfc erfc
 #define gnm_exp exp
diff --git a/src/sf-dpq.c b/src/sf-dpq.c
index 46177e4..6bf9206 100644
--- a/src/sf-dpq.c
+++ b/src/sf-dpq.c
@@ -499,6 +499,18 @@ qlnorm (gnm_float p, gnm_float logmean, gnm_float logsd, gboolean lower_tail, gb
 
 /* ------------------------------------------------------------------------ */
 
+static gnm_float
+dcauchy1 (gnm_float x, const gnm_float shape[], gboolean give_log)
+{
+       return dcauchy (x, shape[0], shape[1], give_log);
+}
+
+static gnm_float
+pcauchy1 (gnm_float x, const gnm_float shape[], gboolean lower_tail, gboolean log_p)
+{
+       return pcauchy (x, shape[0], shape[1], lower_tail, log_p);
+}
+
 /**
  * qcauchy:
  * @p: probability
@@ -515,6 +527,8 @@ gnm_float
 qcauchy (gnm_float p, gnm_float location, gnm_float scale,
         gboolean lower_tail, gboolean log_p)
 {
+       gnm_float x;
+
        if (gnm_isnan(p) || gnm_isnan(location) || gnm_isnan(scale))
                return p + location + scale;
 
@@ -529,13 +543,27 @@ qcauchy (gnm_float p, gnm_float location, gnm_float scale,
                        lower_tail = !lower_tail, p = 0 - gnm_expm1 (p);
                else
                        p = gnm_exp (p);
+               log_p = FALSE;
+       } else {
+               if (p > 0.5) {
+                       p = 1 - p;
+                       lower_tail = !lower_tail;
+               }
        }
-       if (p > 0.5) {
-               p = 1 - p;
-               lower_tail = !lower_tail;
+       x = location + (lower_tail ? -scale : scale) * gnm_cotpi (p);
+
+       if (location != 0 && gnm_abs (x / location) < 0.25) {
+               /* Cancellation has occurred.  */
+               gnm_float shape[2];
+               shape[0] = location;
+               shape[1] = scale;
+               x = pfuncinverter (p, shape, lower_tail, log_p,
+                                  gnm_ninf, gnm_pinf, x,
+                                  pcauchy1, dcauchy1);
+
        }
-       if (lower_tail) scale = -scale;
-       return location + scale / gnm_tanpi (p);
+
+       return x;
 }
 
 /* ------------------------------------------------------------------------ */