[goffice] Quad: fix overflow in mul12.

[Morten Welinder <mortenw src gnome org>]

commit b25fec70148ba52f9a2d886833733f75be67fa36
Author: Morten Welinder <terra gnome org>
Date:   Wed May 17 09:45:02 2017 -0400

    Quad: fix overflow in mul12.
    
    Splitting very large numbers could overflow for no good reason.

 ChangeLog              |    5 +++++
 NEWS                   |    3 +++
 goffice/math/go-quad.c |   35 +++++++++++++++++++++++++++--------
 3 files changed, 35 insertions(+), 8 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 2b0a5fa..2b02e21 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2017-05-17  Morten Welinder  <terra gnome org>
+
+       * goffice/math/go-quad.c (go_quad_mul12): Avoid unnecessary
+       overflow when splitting a number into two parts.
+
 2017-03-20  Morten Welinder <terra gnome org>
 
        * configure.ac: Post-release bump.
diff --git a/NEWS b/NEWS
index 08d4239..61741f6 100644
--- a/NEWS
+++ b/NEWS
@@ -1,5 +1,8 @@
 goffice 0.10.35:
 
+Morten:
+       * Fix quad-precision overflow problem.
+
 --------------------------------------------------------------------------
 goffice 0.10.34:
 
diff --git a/goffice/math/go-quad.c b/goffice/math/go-quad.c
index b018958..9908a86 100644
--- a/goffice/math/go-quad.c
+++ b/goffice/math/go-quad.c
@@ -69,6 +69,7 @@ typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__)));
 #define DOUBLE double
 #define SUFFIX(_n) _n
 #define DOUBLE_MANT_DIG DBL_MANT_DIG
+#define DOUBLE_EPSILON DBL_EPSILON
 
 #ifdef GOFFICE_WITH_LONG_DOUBLE
 #include "go-quad.c"
@@ -76,9 +77,11 @@ typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__)));
 #undef DOUBLE
 #undef SUFFIX
 #undef DOUBLE_MANT_DIG
+#undef DOUBLE_EPSILON
 #define DOUBLE long double
 #define SUFFIX(_n) _n ## l
 #define DOUBLE_MANT_DIG LDBL_MANT_DIG
+#define DOUBLE_EPSILON LDBL_EPSILON
 #endif
 
 #endif
@@ -392,6 +395,22 @@ SUFFIX(go_quad_sub) (QUAD *res, const QUAD *a, const QUAD *b)
        res->l = r - res->h + s;
 }
 
+
+#define SPLIT1(x,h,t) do {                             \
+  DOUBLE p = x * SUFFIX(CST);                          \
+  if (!SUFFIX(go_finite) (p) && SUFFIX(go_finite)(x)) {        \
+     x *= DOUBLE_EPSILON;                              \
+     p = x * SUFFIX(CST);                              \
+     h = x - p + p;                                    \
+     t = x - h;                                                \
+     h *= (1 / DOUBLE_EPSILON);                                \
+     t *= (1 / DOUBLE_EPSILON);                                \
+  } else {                                             \
+     h = x - p + p;                                    \
+     t = x - h;                                                \
+  }                                                    \
+} while (0)
+
 /**
  * go_quad_mul12:
  * @res: (out): result location
@@ -413,20 +432,20 @@ SUFFIX(go_quad_sub) (QUAD *res, const QUAD *a, const QUAD *b)
 void
 SUFFIX(go_quad_mul12) (QUAD *res, DOUBLE x, DOUBLE y)
 {
-       DOUBLE p1 = x * SUFFIX(CST);
-       DOUBLE hx = x - p1 + p1;
-       DOUBLE tx = x - hx;
+       DOUBLE hx, tx, hy, ty, p, q;
 
-       DOUBLE p2 = y * SUFFIX(CST);
-       DOUBLE hy = y - p2 + p2;
-       DOUBLE ty = y - hy;
+       SPLIT1 (x, hx, tx);
+       SPLIT1 (y, hy, ty);
 
-       DOUBLE p = hx * hy;
-       DOUBLE q = hx * ty + tx * hy;
+       p = hx * hy;
+       q = hx * ty + tx * hy;
        res->h = p + q;
        res->l = p - res->h + q + tx * ty;
 }
 
+#undef SPLIT1
+
+
 /**
  * go_quad_mul:
  * @res: (out): result location