Commit d5fb139f authored by Morten Welinder's avatar Morten Welinder

GAMMA: Yet another accuracy tweak.

Also, use go_sinpi etc. from goffice.
parent 9730fe8a
2013-12-19 Morten Welinder <terra@gnome.org>
* src/sf-gamma.c (qfactf): Squeeze a few extra bits out of this,
especially when |x|<1.
2013-12-17 Morten Welinder <terra@gnome.org>
* src/sf-gamma.c (complex_gamma): Turn a complex division into a
......
......@@ -152,7 +152,7 @@ PKG_PROG_PKG_CONFIG(0.18)
dnl *****************************
libspreadsheet_reqs="
libgoffice-${GOFFICE_API_VER} >= 0.10.9
libgoffice-${GOFFICE_API_VER} >= 0.10.10
libgsf-1 >= 1.14.24
libxml-2.0 >= 2.4.12
"
......
......@@ -44,6 +44,7 @@ typedef long double gnm_float;
#define gnm_atanh atanhl
#define gnm_ceil ceill
#define gnm_cosh coshl
#define gnm_cospi go_cospil
#define gnm_erf erfl
#define gnm_erfc erfcl
#define gnm_exp expl
......@@ -75,10 +76,12 @@ typedef long double gnm_float;
#define gnm_pow2 go_pow2l
#define gnm_render_general go_render_generall
#define gnm_sinh sinhl
#define gnm_sinpi go_sinpil
#define gnm_sqrt sqrtl
#define gnm_strto go_strtold
#define gnm_sub_epsilon go_sub_epsilonl
#define gnm_tanh tanhl
#define gnm_tanpi go_tanpil
#ifndef GNM_REDUCES_TRIG_RANGE
#define gnm_sin sinl
#define gnm_cos cosl
......@@ -101,24 +104,30 @@ typedef long double gnm_float;
#define GNM_const(_c) _c ## L
#define GnmQuad GOQuadl
#define gnm_quad_start go_quad_startl
#define gnm_quad_end go_quad_endl
#define gnm_quad_init go_quad_initl
#define gnm_quad_value go_quad_valuel
#define gnm_quad_acos go_quad_acosl
#define gnm_quad_add go_quad_addl
#define gnm_quad_sub go_quad_subl
#define gnm_quad_mul go_quad_mull
#define gnm_quad_asin go_quad_asinl
#define gnm_quad_cos go_quad_cosl
#define gnm_quad_cospi go_quad_cospil
#define gnm_quad_div go_quad_divl
#define gnm_quad_mul12 go_quad_mul12l
#define gnm_quad_sqrt go_quad_sqrtl
#define gnm_quad_pow go_quad_powl
#define gnm_quad_e go_quad_el
#define gnm_quad_end go_quad_endl
#define gnm_quad_exp go_quad_expl
#define gnm_quad_expm1 go_quad_expm1l
#define gnm_quad_zero go_quad_zerol
#define gnm_quad_init go_quad_initl
#define gnm_quad_mul go_quad_mull
#define gnm_quad_mul12 go_quad_mul12l
#define gnm_quad_one go_quad_onel
#define gnm_quad_pi go_quad_pil
#define gnm_quad_e go_quad_el
#define gnm_quad_pow go_quad_powl
#define gnm_quad_sin go_quad_sinl
#define gnm_quad_sinpi go_quad_sinpil
#define gnm_quad_sqrt go_quad_sqrtl
#define gnm_quad_sqrt2 go_quad_sqrt2l
#define gnm_quad_start go_quad_startl
#define gnm_quad_sub go_quad_subl
#define gnm_quad_value go_quad_valuel
#define gnm_quad_zero go_quad_zerol
#define GnmAccumulator GOAccumulatorl
#define gnm_accumulator_start go_accumulator_startl
#define gnm_accumulator_end go_accumulator_endl
......@@ -143,6 +152,7 @@ typedef double gnm_float;
#define gnm_atanh atanh
#define gnm_ceil ceil
#define gnm_cosh cosh
#define gnm_cospi go_cospi
#define gnm_erf erf
#define gnm_erfc erfc
#define gnm_exp exp
......@@ -174,10 +184,12 @@ typedef double gnm_float;
#define gnm_pow2 go_pow2
#define gnm_render_general go_render_general
#define gnm_sinh sinh
#define gnm_sinpi go_sinpi
#define gnm_sqrt sqrt
#define gnm_strto go_strtod
#define gnm_sub_epsilon go_sub_epsilon
#define gnm_tanh tanh
#define gnm_tanpi go_tanpi
#ifndef GNM_REDUCES_TRIG_RANGE
#define gnm_sin sin
#define gnm_cos cos
......@@ -200,24 +212,30 @@ typedef double gnm_float;
#define GNM_const(_c) _c
#define GnmQuad GOQuad
#define gnm_quad_start go_quad_start
#define gnm_quad_end go_quad_end
#define gnm_quad_init go_quad_init
#define gnm_quad_value go_quad_value
#define gnm_quad_acos go_quad_acos
#define gnm_quad_add go_quad_add
#define gnm_quad_sub go_quad_sub
#define gnm_quad_mul go_quad_mul
#define gnm_quad_asin go_quad_asin
#define gnm_quad_cos go_quad_cos
#define gnm_quad_cospi go_quad_cospi
#define gnm_quad_div go_quad_div
#define gnm_quad_mul12 go_quad_mul12
#define gnm_quad_sqrt go_quad_sqrt
#define gnm_quad_pow go_quad_pow
#define gnm_quad_e go_quad_e
#define gnm_quad_end go_quad_end
#define gnm_quad_exp go_quad_exp
#define gnm_quad_expm1 go_quad_expm1
#define gnm_quad_zero go_quad_zero
#define gnm_quad_init go_quad_init
#define gnm_quad_mul go_quad_mul
#define gnm_quad_mul12 go_quad_mul12
#define gnm_quad_one go_quad_one
#define gnm_quad_pi go_quad_pi
#define gnm_quad_e go_quad_e
#define gnm_quad_pow go_quad_pow
#define gnm_quad_sin go_quad_sin
#define gnm_quad_sinpi go_quad_sinpi
#define gnm_quad_sqrt go_quad_sqrt
#define gnm_quad_sqrt2 go_quad_sqrt2
#define gnm_quad_start go_quad_start
#define gnm_quad_sub go_quad_sub
#define gnm_quad_value go_quad_value
#define gnm_quad_zero go_quad_zero
#define GnmAccumulator GOAccumulator
#define gnm_accumulator_start go_accumulator_start
#define gnm_accumulator_end go_accumulator_end
......
......@@ -897,7 +897,8 @@ qfactf (gnm_float x, GnmQuad *mant, int *exp2)
else {
GnmQuad b;
gnm_quad_init (&b, -gnm_sinpi (x)); /* ? */
gnm_quad_init (&b, -x);
gnm_quad_sinpi (&b, &b);
gnm_quad_mul (&b, &b, mant);
gnm_quad_div (mant, &gnm_quad_pi, &b);
*exp2 = -*exp2;
......@@ -946,7 +947,7 @@ qfactf (gnm_float x, GnmQuad *mant, int *exp2)
if (debug) g_printerr ("G(x+1)=%.20g * 2^%d %s\n", gnm_quad_value (mant), *exp2, res ? "overflow" : "");
} else {
GnmQuad s, mFw;
GnmQuad s, qx, mFw;
gnm_float w, f;
int eFw;
......@@ -958,14 +959,13 @@ qfactf (gnm_float x, GnmQuad *mant, int *exp2)
*/
w = gnm_floor (x + 0.5);
f = x - w;
gnm_quad_init (&qx, x);
gnm_quad_init (&s, 1);
while (x < 20) {
GnmQuad a;
x++;
while (w < 20) {
gnm_quad_add (&qx, &qx, &gnm_quad_one);
w++;
gnm_quad_init (&a, x);
gnm_quad_mul (&s, &s, &a);
gnm_quad_mul (&s, &s, &qx);
}
if (qfacti ((int)w, &mFw, &eFw)) {
......
......@@ -78,88 +78,6 @@ gnm_acoth (gnm_float x)
/* ------------------------------------------------------------------------- */
/**
* gnm_sinpi:
* @x: a number
*
* Returns: the sine of Pi times @x, but with less error than doing the
* multiplication outright.
*/
gnm_float
gnm_sinpi (gnm_float x)
{
int k;
if (gnm_isnan (x))
return x;
if (!gnm_finite (x))
return gnm_nan;
k = (x < 0) ? 2 : 0;
x = gnm_fmod (gnm_abs (x), 2);
if (x >= 1) {
x -= 1;
k ^= 2;
}
if (x >= 0.5) {
x -= 0.5;
k += 1;
}
if (x == 0) {
static const gnm_float ys[4] = { 0, 1, -0, -1 };
return ys[k];
} else {
switch (k) {
default:
case 0: return +gnm_sin (M_PIgnum * x);
case 1: return +gnm_cos (M_PIgnum * x);
case 2: return -gnm_sin (M_PIgnum * x);
case 3: return -gnm_cos (M_PIgnum * x);
}
}
}
/**
* gnm_cospi:
* @x: a number
*
* Returns: the cosine of Pi times @x, but with less error than doing the
* multiplication outright.
*/
gnm_float
gnm_cospi (gnm_float x)
{
int k = 0;
if (!gnm_finite (x))
return gnm_nan;
x = gnm_fmod (gnm_abs (x), 2);
if (x >= 1) {
x -= 1;
k ^= 2;
}
if (x >= 0.5) {
x -= 0.5;
k += 1;
}
if (x == 0) {
static const gnm_float ys[4] = { 1, 0, -1, -0 };
return ys[k];
} else {
switch (k) {
default:
case 0: return +gnm_cos (M_PIgnum * x);
case 1: return -gnm_sin (M_PIgnum * x);
case 2: return -gnm_cos (M_PIgnum * x);
case 3: return +gnm_sin (M_PIgnum * x);
}
}
}
/* ------------------------------------------------------------------------- */
#ifdef GNM_REDUCES_TRIG_RANGE
......
......@@ -8,9 +8,6 @@ gnm_float gnm_acot (gnm_float x);
gnm_float gnm_coth (gnm_float x);
gnm_float gnm_acoth (gnm_float x);
gnm_float gnm_sinpi (gnm_float x);
gnm_float gnm_cospi (gnm_float x);
#ifdef GNM_REDUCES_TRIG_RANGE
/* gnm_sin, gnm_cos, gnm_tan prototyped in numbers.h */
#endif
......
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