[goffice] Quad: add hypot.

[Morten Welinder <mortenw src gnome org>]

commit 9863b3d57c79c7be40155d2aa38ff6816cc03f98
Author: Morten Welinder <terra gnome org>
Date:   Mon Dec 9 14:30:16 2013 -0500

    Quad: add hypot.

 NEWS                      |    2 +-
 goffice/math/go-complex.c |    5 +----
 goffice/math/go-quad.c    |   34 +++++++++++++++++++++++++++++++++-
 goffice/math/go-quad.h    |    2 ++
 4 files changed, 37 insertions(+), 6 deletions(-)
---
diff --git a/NEWS b/NEWS
index c21d8d4..2944e00 100644
--- a/NEWS
+++ b/NEWS
@@ -5,7 +5,7 @@ Jean
 
 Morten:
        * Improve complex math, notably complex power.
-       * Add quad precision log, atan2 functions.
+       * Add quad precision log, atan2, hypot functions.
 
 --------------------------------------------------------------------------
 goffice 0.10.9:
diff --git a/goffice/math/go-complex.c b/goffice/math/go-complex.c
index 24b3900..ceb3501 100644
--- a/goffice/math/go-complex.c
+++ b/goffice/math/go-complex.c
@@ -388,10 +388,7 @@ SUFFIX(go_complex_pow) (COMPLEX *dst, COMPLEX const *a, COMPLEX const *b)
                SUFFIX(go_quad_init) (&qa, a->re);
                SUFFIX(go_quad_init) (&qb, a->im);
                SUFFIX(go_quad_atan2pi) (&qarg, &qb, &qa);
-               SUFFIX(go_quad_mul) (&qa, &qa, &qa);
-               SUFFIX(go_quad_mul) (&qb, &qb, &qb);
-               SUFFIX(go_quad_add) (&qa, &qa, &qb);
-               SUFFIX(go_quad_sqrt) (&qr, &qa);
+               SUFFIX(go_quad_hypot) (&qr, &qa, &qb);
 
                /*
                 * This is the square root of the power we really want,
diff --git a/goffice/math/go-quad.c b/goffice/math/go-quad.c
index d63b529..2b2f114 100644
--- a/goffice/math/go-quad.c
+++ b/goffice/math/go-quad.c
@@ -618,7 +618,7 @@ SUFFIX(go_quad_pow_frac) (QUAD *res, const QUAD *x, const QUAD *y,
 /*
  * This computes pow(x,y) in the following way:
  *
- * 1. y is considered the sum of a number of one-but values.
+ * 1. y is considered the sum of a number of one-bit values.
  * 2. For each y bit in the integer part of y, the corresponding x^y
  *    is computed by repeated squaring.
  * 3. For each y bit in the fractional part of y, the corresponding x^y
@@ -922,3 +922,35 @@ SUFFIX(go_quad_atan2pi) (QUAD *res, const QUAD *y, const QUAD *x)
        SUFFIX(go_quad_atan2) (res, y, x);
        SUFFIX(go_quad_div) (res, res, &SUFFIX(go_quad_pi));
 }
+
+void
+SUFFIX(go_quad_hypot) (QUAD *res, const QUAD *a, const QUAD *b)
+{
+       int e;
+       QUAD qa2, qb2, qn;
+
+       if (a->h == 0) {
+               *res = *b;
+               return;
+       }
+       if (b->h == 0) {
+               *res = *a;
+               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);
+}
diff --git a/goffice/math/go-quad.h b/goffice/math/go-quad.h
index 912b76e..e0d7cf5 100644
--- a/goffice/math/go-quad.h
+++ b/goffice/math/go-quad.h
@@ -29,6 +29,7 @@ void go_quad_expm1 (GOQuad *res, const GOQuad *a);
 void go_quad_log (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_hypot (GOQuad *res, const GOQuad *a, const GOQuad *b);
 
 void go_quad_mul12 (GOQuad *res, double x, double y);
 
@@ -70,6 +71,7 @@ void go_quad_expm1l (GOQuadl *res, const GOQuadl *a);
 void go_quad_logl (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_hypotl (GOQuadl *res, const GOQuadl *a, const GOQuadl *b);
 
 void go_quad_mul12l (GOQuadl *res, long double x, long double y);