[goffice] Quad: add sin, cos.

[Morten Welinder <mortenw src gnome org>]

commit 3efd441af713625fe8936bd947e9c3b85eca1b55
Author: Morten Welinder <terra gnome org>
Date:   Wed Dec 18 19:13:31 2013 -0500

    Quad: add sin, cos.

 NEWS                   |    3 +-
 goffice/math/go-quad.c |  147 +++++++++++++++++++++++++++++++++++++++++++++--
 goffice/math/go-quad.h |   12 +++-
 tests/test-quad.c      |   12 ++++
 4 files changed, 162 insertions(+), 12 deletions(-)
---
diff --git a/NEWS b/NEWS
index 65f3726..44284a6 100644
--- a/NEWS
+++ b/NEWS
@@ -7,7 +7,8 @@ Jean
 
 Morten:
        * Improve complex math, notably complex power.
-       * Add quad precision log, atan2, hypot, asin, acos functions.
+       * Add quad precision log and hypot functions.
+       * Add quad precision atan2, sin, asin, cos, acos functions.
        * Add go_sinpi, go_cospi, go_tanpi, go_atan2pi.
 
 --------------------------------------------------------------------------
diff --git a/goffice/math/go-quad.c b/goffice/math/go-quad.c
index 83cb823..af94e38 100644
--- a/goffice/math/go-quad.c
+++ b/goffice/math/go-quad.c
@@ -807,6 +807,18 @@ SUFFIX(go_quad_hypot) (QUAD *res, const QUAD *a, const QUAD *b)
        res->l = SUFFIX(ldexp) (qn.l, e);
 }
 
+/* sqrt(1-a*a) helper */
+static void
+SUFFIX(go_quad_ihypot) (QUAD *res, const QUAD *a)
+{
+       QUAD qp;
+
+       SUFFIX(go_quad_mul) (&qp, a, a);
+       SUFFIX(go_quad_sub) (&qp, &SUFFIX(go_quad_one), &qp);
+       SUFFIX(go_quad_sqrt) (res, &qp);
+}
+
+
 #ifdef DEFINE_COMMON
 typedef enum {
        AGM_ARCSIN,
@@ -836,9 +848,7 @@ SUFFIX(go_quad_agm_internal) (QUAD *res, AGM_Method method, const QUAD *x)
        switch (method) {
        case AGM_ARCSIN:
                g_return_if_fail (SUFFIX(fabs) (x->h) <= 1);
-               SUFFIX(go_quad_mul) (&dpk[0], x, x);
-               SUFFIX(go_quad_sub) (&dpk[0], &SUFFIX(go_quad_one), &dpk[0]);
-               SUFFIX(go_quad_sqrt) (&dpk[0], &dpk[0]);
+               SUFFIX(go_quad_ihypot) (&dpk[0], x);
                gp = SUFFIX(go_quad_one);
                qnum = *x;
                break;
@@ -846,9 +856,7 @@ SUFFIX(go_quad_agm_internal) (QUAD *res, AGM_Method method, const QUAD *x)
                g_return_if_fail (SUFFIX(fabs) (x->h) <= 1);
                dpk[0] = *x;
                gp = SUFFIX(go_quad_one);
-               SUFFIX(go_quad_mul) (&qnum, x, x);
-               SUFFIX(go_quad_sub) (&qnum, &SUFFIX(go_quad_one), &qnum);
-               SUFFIX(go_quad_sqrt) (&qnum, &qnum);
+               SUFFIX(go_quad_ihypot) (&qnum, x);
                break;
        case AGM_ARCTAN:
                g_return_if_fail (SUFFIX(fabs) (x->h) <= 1);
@@ -879,7 +887,7 @@ SUFFIX(go_quad_agm_internal) (QUAD *res, AGM_Method method, const QUAD *x)
 
                SUFFIX(go_quad_div) (&qr, &qnum, &dk[n]);
                SUFFIX(go_quad_sub) (&qrp, &qrp, &qr);
-               if (SUFFIX(fabs)(qrp.h) <= SUFFIX(ldexp) (SUFFIX(fabs)(qr.h), -2 * (DBL_MANT_DIG - 1))) {
+               if (SUFFIX(fabs)(qrp.h) <= SUFFIX(ldexp) (SUFFIX(fabs)(qr.h), -2 * (DOUBLE_MANT_DIG - 1))) {
                        converged = TRUE;
                        break;
                }
@@ -984,6 +992,113 @@ SUFFIX(go_quad_atan2pi) (QUAD *res, const QUAD *y, const QUAD *x)
        SUFFIX(go_quad_div) (res, res, &SUFFIX(go_quad_pi));
 }
 
+static gboolean
+SUFFIX(reduce_pi_half) (QUAD *res, const QUAD *a, int *pk)
+{
+       static QUAD pi_half;
+       QUAD qa, qk, qh, qb;
+       DOUBLE k;
+       unsigned ui;
+       static const DOUBLE pi_half_parts[] = {
+               +0x1.921fb54442d18p+0,
+               +0x1.1a62633145c04p-54,
+               +0x1.707344a40938p-105,
+               +0x1.114cf98e80414p-156,
+               +0x1.bea63b139b224p-207,
+               +0x1.14a08798e3404p-259,
+               +0x1.bbdf2a33679a4p-312,
+               +0x1.a431b302b0a6cp-363,
+               +0x1.f25f14374fe1p-415,
+               +0x1.ab6b6a8e122fp-466
+       };
+
+       if (!SUFFIX(go_finite) (a->h))
+               return TRUE;
+
+       if (SUFFIX(fabs) (a->h) > SUFFIX(ldexp) (1.0, DOUBLE_MANT_DIG)) {
+               g_warning ("Reduced accuracy for very large trigonometric arguments");
+               return TRUE;
+       }
+
+       if (pi_half.h == 0) {
+               pi_half.h = SUFFIX(go_quad_pi).h * 0.5;
+               pi_half.l = SUFFIX(go_quad_pi).l * 0.5;
+       }
+
+       SUFFIX(go_quad_div) (&qk, a, &pi_half);
+       qh.h = 0.5; qh.l = 0;
+       SUFFIX(go_quad_add) (&qk, &qk, &qh);
+       SUFFIX(go_quad_floor) (&qk, &qk);
+       k = SUFFIX(go_quad_value) (&qk);
+       *pk = (int)(SUFFIX(fmod) (k, 4));
+
+       qa = *a;
+       for (ui = 0; ui < G_N_ELEMENTS(pi_half_parts); ui++) {
+               SUFFIX(go_quad_mul12) (&qb, pi_half_parts[ui], k);
+               SUFFIX(go_quad_sub) (&qa, &qa, &qb);
+       }
+
+       *res = qa;
+
+       return FALSE;
+}
+
+static void
+SUFFIX(do_sin) (QUAD *res, const QUAD *a, int k)
+{
+       QUAD qr;
+
+       if (k & 1) {
+               QUAD qn, qd, qq, aa;
+
+               aa.h = SUFFIX(fabs)(a->h);
+               aa.l = SUFFIX(fabs)(a->l);
+               SUFFIX(go_quad_init) (&qr, cos (aa.h));
+
+               /* Newton step */
+               SUFFIX(go_quad_acos) (&qn, &qr);
+               SUFFIX(go_quad_sub) (&qn, &qn, &aa);
+               SUFFIX(go_quad_ihypot) (&qd, &qr);
+               SUFFIX(go_quad_mul) (&qq, &qn, &qd);
+               SUFFIX(go_quad_add) (&qr, &qr, &qq);
+       } else {
+               QUAD qn, qd, qq;
+               SUFFIX(go_quad_init) (&qr, sin (a->h));
+
+               /* Newton step */
+               SUFFIX(go_quad_asin) (&qn, &qr);
+               SUFFIX(go_quad_sub) (&qn, &qn, a);
+               SUFFIX(go_quad_ihypot) (&qd, &qr);
+               SUFFIX(go_quad_mul) (&qq, &qn, &qd);
+               SUFFIX(go_quad_sub) (&qr, &qr, &qq);
+       }
+
+       if (k & 2) {
+               qr.h = 0 - qr.h;
+               qr.l = 0 - qr.l;
+       }
+
+       *res = qr;
+}
+
+/**
+ * go_quad_sin: (skip)
+ **/
+/**
+ * go_quad_sinl: (skip)
+ **/
+void
+SUFFIX(go_quad_sin) (QUAD *res, const QUAD *a)
+{
+       int k;
+       QUAD a0;
+
+       if (SUFFIX(reduce_pi_half) (&a0, a, &k))
+               SUFFIX(go_quad_init) (res, SUFFIX(sin) (a->h));
+       else
+               SUFFIX(do_sin) (res, &a0, k);
+}
+
 /**
  * go_quad_asin: (skip)
  **/
@@ -1007,6 +1122,24 @@ SUFFIX(go_quad_asin) (QUAD *res, const QUAD *a)
 }
 
 /**
+ * go_quad_cos: (skip)
+ **/
+/**
+ * go_quad_cosl: (skip)
+ **/
+void
+SUFFIX(go_quad_cos) (QUAD *res, const QUAD *a)
+{
+       int k;
+       QUAD a0;
+
+       if (SUFFIX(reduce_pi_half) (&a0, a, &k))
+               SUFFIX(go_quad_init) (res, SUFFIX(cos) (a->h));
+       else
+               SUFFIX(do_sin) (res, &a0, k + 1);
+}
+
+/**
  * go_quad_acos: (skip)
  **/
 /**
diff --git a/goffice/math/go-quad.h b/goffice/math/go-quad.h
index 8227b2a..99a2607 100644
--- a/goffice/math/go-quad.h
+++ b/goffice/math/go-quad.h
@@ -29,10 +29,12 @@ void go_quad_expm1 (GOQuad *res, const GOQuad *a);
 void go_quad_log (GOQuad *res, const GOQuad *a);
 void go_quad_hypot (GOQuad *res, const GOQuad *a, const GOQuad *b);
 
-void go_quad_atan2 (GOQuad *res, const GOQuad *y, const GOQuad *x);
-void go_quad_atan2pi (GOQuad *res, const GOQuad *y, const GOQuad *x);
+void go_quad_sin (GOQuad *res, const GOQuad *a);
 void go_quad_asin (GOQuad *res, const GOQuad *a);
+void go_quad_cos (GOQuad *res, const GOQuad *a);
 void go_quad_acos (GOQuad *res, const GOQuad *a);
+void go_quad_atan2 (GOQuad *res, const GOQuad *y, const GOQuad *x);
+void go_quad_atan2pi (GOQuad *res, const GOQuad *y, const GOQuad *x);
 
 void go_quad_mul12 (GOQuad *res, double x, double y);
 
@@ -74,10 +76,12 @@ void go_quad_expm1l (GOQuadl *res, const GOQuadl *a);
 void go_quad_logl (GOQuadl *res, const GOQuadl *a);
 void go_quad_hypotl (GOQuadl *res, const GOQuadl *a, const GOQuadl *b);
 
-void go_quad_atan2l (GOQuadl *res, const GOQuadl *y, const GOQuadl *x);
-void go_quad_atan2pil (GOQuadl *res, const GOQuadl *y, const GOQuadl *x);
+void go_quad_sinl (GOQuadl *res, const GOQuadl *a);
 void go_quad_asinl (GOQuadl *res, const GOQuadl *a);
+void go_quad_cosl (GOQuadl *res, const GOQuadl *a);
 void go_quad_acosl (GOQuadl *res, const GOQuadl *a);
+void go_quad_atan2l (GOQuadl *res, const GOQuadl *y, const GOQuadl *x);
+void go_quad_atan2pil (GOQuadl *res, const GOQuadl *y, const GOQuadl *x);
 
 void go_quad_mul12l (GOQuadl *res, long double x, long double y);
 
diff --git a/tests/test-quad.c b/tests/test-quad.c
index 88f8514..beed1b3 100644
--- a/tests/test-quad.c
+++ b/tests/test-quad.c
@@ -289,9 +289,16 @@ floor_tests (void)
        UNTEST1(a_,go_quad_acos,acos,"acos");   \
 } while (0)
 
+#define TEST2(a_) do {                         \
+       UNTEST1(a_,go_quad_sin,sin,"sin");      \
+       UNTEST1(a_,go_quad_cos,cos,"cos");      \
+} while (0)
+
 static void
 trig_tests (void)
 {
+       double d;
+
        TEST1 (0);
        TEST1 (0.25);
        TEST1 (0.5);
@@ -301,6 +308,11 @@ trig_tests (void)
        TEST1 (-0.5);
        TEST1 (-0.75);
        TEST1 (-1);
+
+       for (d = 0; d < 10; d += 0.125) {
+               TEST2(d);
+               TEST2(-d);
+       }
 }
 
 #undef TEST1