[gnumeric] GAMMA: Yet another accuracy tweak.

[Morten Welinder <mortenw src gnome org>]

commit d5fb139f0e11776951ebad0451dcb88a7a6396f3
Author: Morten Welinder <terra gnome org>
Date:   Thu Dec 19 13:14:02 2013 -0500

    GAMMA: Yet another accuracy tweak.
    
    Also, use go_sinpi etc. from goffice.

 ChangeLog      |    5 +++
 configure.ac   |    2 +-
 src/numbers.h  |   62 +++++++++++++++++++++++++++---------------
 src/sf-gamma.c |   14 +++++-----
 src/sf-trig.c  |   82 --------------------------------------------------------
 src/sf-trig.h  |    3 --
 6 files changed, 53 insertions(+), 115 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index be319a7..6f4afd9 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+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
diff --git a/configure.ac b/configure.ac
index 04048d0..46c7d8b 100644
--- a/configure.ac
+++ b/configure.ac
@@ -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
 "
diff --git a/src/numbers.h b/src/numbers.h
index 3a74cb4..cc09273 100644
--- a/src/numbers.h
+++ b/src/numbers.h
@@ -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
diff --git a/src/sf-gamma.c b/src/sf-gamma.c
index f03510e..66030c0 100644
--- a/src/sf-gamma.c
+++ b/src/sf-gamma.c
@@ -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)) {
diff --git a/src/sf-trig.c b/src/sf-trig.c
index 93e3d1b..3ad232a 100644
--- a/src/sf-trig.c
+++ b/src/sf-trig.c
@@ -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
 
diff --git a/src/sf-trig.h b/src/sf-trig.h
index c9ad770..fb1c72c 100644
--- a/src/sf-trig.h
+++ b/src/sf-trig.h
@@ -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