Commit 42e56b0a authored by Morten Welinder's avatar Morten Welinder

Math: improve accuracy of dnorm.

Doing exp(-x*x/2) is not accurate for large x.  Large starts about 5.
This affects a large number of functions that use dnorm.
parent 83097bcc
2013-12-30 Morten Welinder <terra@gnome.org>
* src/mathfunc.c (dnorm): Improve accuracy for x>5 (normalized).
2013-12-25 Morten Welinder <terra@gnome.org>
* src/item-grid.c (cb_cursor_come_to_rest): Clear ->tip_timer when
......
......@@ -21,6 +21,7 @@ Morten:
* Fix timeout related critical in ItemGrid.
* Fix win32 print crash. [#719997]
* Fix solver problem with constraints.
* Improve accuracy of r.PNORM for large arguments.
--------------------------------------------------------------------------
Gnumeric 1.12.9
......
......@@ -234,12 +234,26 @@ gnm_float dnorm(gnm_float x, gnm_float mu, gnm_float sigma, gboolean give_log)
return (x == mu) ? gnm_pinf : R_D__0;
}
x = (x - mu) / sigma;
if(!gnm_finite(x)) return R_D__0;
return (give_log ?
-(M_LN_SQRT_2PI + 0.5 * x * x + gnm_log(sigma)) :
M_1_SQRT_2PI * gnm_exp(-0.5 * x * x) / sigma);
/* M_1_SQRT_2PI = 1 / gnm_sqrt(2 * pi) */
x = gnm_abs (x);
if (give_log)
return -(M_LN_SQRT_2PI + 0.5 * x * x + gnm_log(sigma));
else if (x < 5)
return M_1_SQRT_2PI * gnm_exp(-0.5 * x * x) / sigma;
else if (x >= 256)
return 0; /* Will underflow anyway. */
else {
/*
* Split x into two parts, x=x1+x2, such that x2 is
* small and x1 has less than 26 bits. That ensures
* that x1*x1 is error free.
*/
gnm_float x1 = gnm_floor (x * 65536 + 0.5) / 65536;
gnm_float x2 = x - x1;
return M_1_SQRT_2PI / sigma *
(gnm_exp(-0.5 * x1 * x1) *
gnm_exp(-(0.5 * x2 + x1) * x2));
}
}
/* ------------------------------------------------------------------------ */
......
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