Commit a838504d authored by Morten Welinder's avatar Morten Welinder

R.QCAUCHY: Improve accuracy.

With non-zero location we can suffer cancellation.
parent fbbff099
......@@ -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.
......
......@@ -33,6 +33,7 @@ Morten:
* Fix sheet filter problem.
* Minor R.DNORM improvement.
* Improve RAYLEIGH accuracy.
* Improve R.QCAUCHY accuracy.
--------------------------------------------------------------------------
Gnumeric 1.12.21
......
......@@ -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));
}
/* ------------------------------------------------------------------------ */
......
......@@ -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
......
......@@ -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;
}
/* ------------------------------------------------------------------------ */
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment