[goffice] Quad: add asin and acos.

[Morten Welinder <mortenw src gnome org>]

commit a5980b264e458c6251fab01bfad8b42426c9bc11
Author: Morten Welinder <terra gnome org>
Date:   Tue Dec 17 14:59:41 2013 -0500

    Quad: add asin and acos.

 ChangeLog              |    5 ++
 NEWS                   |    2 +-
 goffice/math/go-quad.c |  149 +++++++++++++++++++++++++++++++++++------------
 goffice/math/go-quad.h |   10 +++-
 tests/test-quad.c      |   26 ++++++++
 5 files changed, 151 insertions(+), 41 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 654ab94..6758a7c 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,10 @@
 2013-12-17  Morten Welinder  <terra gnome org>
 
+       * tests/test-quad.c (floor_tests): Test asin and acos.
+
+       * goffice/math/go-quad.c (go_quad_acos, go_quad_asin): new
+       functions.
+
        * goffice/math/go-complex.c (go_complex_from_polar_pi)
        (go_complex_angle_pi): Handle diagonals too.
 
diff --git a/NEWS b/NEWS
index beab22e..eb1cd54 100644
--- a/NEWS
+++ b/NEWS
@@ -7,7 +7,7 @@ Jean
 
 Morten:
        * Improve complex math, notably complex power.
-       * Add quad precision log, atan2, hypot functions.
+       * Add quad precision log, atan2, hypot, asin, acos functions.
 
 --------------------------------------------------------------------------
 goffice 0.10.9:
diff --git a/goffice/math/go-quad.c b/goffice/math/go-quad.c
index 8e9c520..83cb823 100644
--- a/goffice/math/go-quad.c
+++ b/goffice/math/go-quad.c
@@ -773,30 +773,92 @@ SUFFIX(go_quad_log) (QUAD *res, const QUAD *a)
        }
 }
 
+void
+SUFFIX(go_quad_hypot) (QUAD *res, const QUAD *a, const QUAD *b)
+{
+       int e;
+       QUAD qa2, qb2, qn;
+
+       if (a->h == 0) {
+               res->h = SUFFIX(fabs)(b->h);
+               res->l = SUFFIX(fabs)(b->l);
+               return;
+       }
+       if (b->h == 0) {
+               res->h = SUFFIX(fabs)(a->h);
+               res->l = SUFFIX(fabs)(a->l);
+               return;
+       }
+
+       /* Scale by power of 2 to protect against over- and underflow */
+       (void)SUFFIX(frexp) (MAX (SUFFIX(fabs) (a->h), SUFFIX(fabs) (b->h)), &e);
+
+       qa2.h = SUFFIX(ldexp) (a->h, -e);
+       qa2.l = SUFFIX(ldexp) (a->l, -e);
+       SUFFIX(go_quad_mul) (&qa2, &qa2, &qa2);
+
+       qb2.h = SUFFIX(ldexp) (b->h, -e);
+       qb2.l = SUFFIX(ldexp) (b->l, -e);
+       SUFFIX(go_quad_mul) (&qb2, &qb2, &qb2);
+
+       SUFFIX(go_quad_add) (&qn, &qa2, &qb2);
+       SUFFIX(go_quad_sqrt) (&qn, &qn);
+       res->h = SUFFIX(ldexp) (qn.h, e);
+       res->l = SUFFIX(ldexp) (qn.l, e);
+}
+
+#ifdef DEFINE_COMMON
+typedef enum {
+       AGM_ARCSIN,
+       AGM_ARCCOS,
+       AGM_ARCTAN
+} AGM_Method;
+#endif
+
 static void
-SUFFIX(go_quad_atan_internal) (QUAD *res, const QUAD *x)
+SUFFIX(go_quad_agm_internal) (QUAD *res, AGM_Method method, const QUAD *x)
 {
-       QUAD g, gp, dk[20], dpk[20], qr, qrp;
+       QUAD g, gp, dk[20], dpk[20], qr, qrp, qnum;
        int n, k;
        gboolean converged = FALSE;
 
-       g_return_if_fail (SUFFIX(fabs) (x->h) <= 1);
-
        /*
         * This follows "An Algorithm for Computing Logarithms
         * and Arctangents" by B. C. Carlson in *Mathematics of
         * Computation*, Volume 26, Number 118, April 1972.
         *
-        * If need be we can do log, arcsin, arccos, arctanh,
-        * arcsinh, and arccosh we the same code.
+        * If need be we can do log, arctanh, arcsinh, and
+        * arccosh we the same code.
         */
 
        qrp = SUFFIX(go_quad_zero);
 
-       dpk[0] = SUFFIX(go_quad_one);
-       SUFFIX(go_quad_mul) (&gp, x, x);
-       SUFFIX(go_quad_add) (&gp, &gp, &SUFFIX(go_quad_one));
-       SUFFIX(go_quad_sqrt) (&gp, &gp);
+       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]);
+               gp = SUFFIX(go_quad_one);
+               qnum = *x;
+               break;
+       case AGM_ARCCOS:
+               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);
+               break;
+       case AGM_ARCTAN:
+               g_return_if_fail (SUFFIX(fabs) (x->h) <= 1);
+               dpk[0] = SUFFIX(go_quad_one);
+               SUFFIX(go_quad_hypot) (&gp, x, &SUFFIX(go_quad_one));
+               qnum = *x;
+               break;
+       default:
+               g_assert_not_reached ();
+       }
 
        for (n = 1; n < (int)G_N_ELEMENTS(dk); n++) {
                SUFFIX(go_quad_add) (&dk[0], &dpk[0], &gp);
@@ -815,8 +877,7 @@ SUFFIX(go_quad_atan_internal) (QUAD *res, const QUAD *x)
                        SUFFIX(go_quad_div) (&dk[k], &dk[k], &f);
                }
 
-               SUFFIX(go_quad_div) (&qr, x, &dk[n]);
-               
+               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))) {
                        converged = TRUE;
@@ -829,7 +890,7 @@ SUFFIX(go_quad_atan_internal) (QUAD *res, const QUAD *x)
        }
 
        if (!converged)
-               g_warning ("go_quad_atan_internal(%.20g) failed to converge\n",
+               g_warning ("go_quad_agm_internal(%.20g) failed to converge\n",
                           (double)SUFFIX(go_quad_value) (x));
 
        *res = qr;
@@ -879,13 +940,13 @@ SUFFIX(go_quad_atan2) (QUAD *res, const QUAD *y, const QUAD *x)
 
        if (SUFFIX(fabs) (dx) >= SUFFIX(fabs) (dy)) {
                SUFFIX(go_quad_div) (&qr, y, x);
-               SUFFIX(go_quad_atan_internal) (res, &qr);
+               SUFFIX(go_quad_agm_internal) (res, AGM_ARCTAN, &qr);
        } else {
                DOUBLE f;
                QUAD qa;
 
                SUFFIX(go_quad_div) (&qr, x, y);
-               SUFFIX(go_quad_atan_internal) (res, &qr);
+               SUFFIX(go_quad_agm_internal) (res, AGM_ARCTAN, &qr);
 
                f = (qr.h >= 0) ? 0.5 : -0.5;
                qa.h = f * SUFFIX(go_quad_pi).h;
@@ -923,36 +984,48 @@ SUFFIX(go_quad_atan2pi) (QUAD *res, const QUAD *y, const QUAD *x)
        SUFFIX(go_quad_div) (res, res, &SUFFIX(go_quad_pi));
 }
 
+/**
+ * go_quad_asin: (skip)
+ **/
+/**
+ * go_quad_asinl: (skip)
+ **/
 void
-SUFFIX(go_quad_hypot) (QUAD *res, const QUAD *a, const QUAD *b)
+SUFFIX(go_quad_asin) (QUAD *res, const QUAD *a)
 {
-       int e;
-       QUAD qa2, qb2, qn;
+       QUAD aa, aam1;
 
-       if (a->h == 0) {
-               res->h = SUFFIX(fabs)(b->h);
-               res->l = SUFFIX(fabs)(b->l);
-               return;
-       }
-       if (b->h == 0) {
-               res->h = SUFFIX(fabs)(a->h);
-               res->l = SUFFIX(fabs)(a->l);
+       aa.h = SUFFIX(fabs) (a->h);
+       aa.l = SUFFIX(fabs) (a->l);
+       SUFFIX(go_quad_sub) (&aam1, &aa, &SUFFIX(go_quad_one));
+       if (aam1.h > 0) {
+               SUFFIX(go_quad_init) (res, SUFFIX(go_nan));
                return;
        }
 
-       /* Scale by power of 2 to protect against over- and underflow */
-       (void)SUFFIX(frexp) (MAX (SUFFIX(fabs) (a->h), SUFFIX(fabs) (b->h)), &e);
+       SUFFIX(go_quad_agm_internal) (res, AGM_ARCSIN, a);
+}
 
-       qa2.h = SUFFIX(ldexp) (a->h, -e);
-       qa2.l = SUFFIX(ldexp) (a->l, -e);
-       SUFFIX(go_quad_mul) (&qa2, &qa2, &qa2);
+/**
+ * go_quad_acos: (skip)
+ **/
+/**
+ * go_quad_acosl: (skip)
+ **/
+void
+SUFFIX(go_quad_acos) (QUAD *res, const QUAD *a)
+{
+       QUAD aa, aam1;
 
-       qb2.h = SUFFIX(ldexp) (b->h, -e);
-       qb2.l = SUFFIX(ldexp) (b->l, -e);
-       SUFFIX(go_quad_mul) (&qb2, &qb2, &qb2);
+       aa.h = SUFFIX(fabs) (a->h);
+       aa.l = SUFFIX(fabs) (a->l);
+       SUFFIX(go_quad_sub) (&aam1, &aa, &SUFFIX(go_quad_one));
+       if (aam1.h > 0) {
+               SUFFIX(go_quad_init) (res, SUFFIX(go_nan));
+               return;
+       }
 
-       SUFFIX(go_quad_add) (&qn, &qa2, &qb2);
-       SUFFIX(go_quad_sqrt) (&qn, &qn);
-       res->h = SUFFIX(ldexp) (qn.h, e);
-       res->l = SUFFIX(ldexp) (qn.l, e);
+       SUFFIX(go_quad_agm_internal) (res, AGM_ARCCOS, &aa);
+       if (a->h < 0)
+               SUFFIX(go_quad_sub) (res, &SUFFIX(go_quad_pi), res);
 }
diff --git a/goffice/math/go-quad.h b/goffice/math/go-quad.h
index e0d7cf5..8227b2a 100644
--- a/goffice/math/go-quad.h
+++ b/goffice/math/go-quad.h
@@ -27,9 +27,12 @@ void go_quad_pow (GOQuad *res, double *exp2, const GOQuad *x, const GOQuad *y);
 void go_quad_exp (GOQuad *res, double *exp2, const GOQuad *a);
 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_hypot (GOQuad *res, const GOQuad *a, const GOQuad *b);
+void go_quad_asin (GOQuad *res, const GOQuad *a);
+void go_quad_acos (GOQuad *res, const GOQuad *a);
 
 void go_quad_mul12 (GOQuad *res, double x, double y);
 
@@ -69,9 +72,12 @@ void go_quad_powl (GOQuadl *res, long double *exp2, const GOQuadl *x, const GOQu
 void go_quad_expl (GOQuadl *res, long double *exp2, const GOQuadl *a);
 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_hypotl (GOQuadl *res, const GOQuadl *a, const GOQuadl *b);
+void go_quad_asinl (GOQuadl *res, const GOQuadl *a);
+void go_quad_acosl (GOQuadl *res, const GOQuadl *a);
 
 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 e564b09..b4fe2f9 100644
--- a/tests/test-quad.c
+++ b/tests/test-quad.c
@@ -277,6 +277,31 @@ floor_tests (void)
        g_assert (go_quad_value (&b) == 0.5);
 }
 
+#undef TEST1
+
+/* ------------------------------------------------------------------------- */
+
+#define TEST1(a_) do {                         \
+       UNTEST1(a_,go_quad_asin,asin,"asin");   \
+       UNTEST1(a_,go_quad_acos,acos,"acos");   \
+} while (0)
+
+static void
+trig_tests (void)
+{
+       TEST1 (0);
+       TEST1 (0.25);
+       TEST1 (0.5);
+       TEST1 (0.75);
+       TEST1 (1);
+       TEST1 (-0.25);
+       TEST1 (-0.5);
+       TEST1 (-0.75);
+       TEST1 (-1);
+}
+
+#undef TEST1
+
 /* ------------------------------------------------------------------------- */
 
 static void
@@ -327,6 +352,7 @@ main (int argc, char **argv)
        floor_tests ();
        atan2_tests ();
        hypot_tests ();
+       trig_tests ();
 
        go_quad_end (state);