This is the mail archive of the
gcc-patches@gcc.gnu.org
mailing list for the GCC project.
[PATCH] Precision tracking in real.c
- From: Roger Sayle <roger at www dot eyesopen dot com>
- To: gcc-patches at gcc dot gnu dot org
- Date: Sun, 20 Apr 2003 14:09:59 -0600 (MDT)
- Subject: [PATCH] Precision tracking in real.c
Several maintainers have voiced concern about the accuracy of GCC's
floating point optimizations. The following patch makes it to possible
to track whether a floating point result is guaranteed to be exact, for
example that pow(3.0,4.0) is identically equal to 81.0. This bookkeeping
allows such constant folding transformations to be applied even without
flag_unsafe_math_optimizations as the result is guaranteed to be atleast
as accurate as the system math library, which often does much worse.
There are no functional changes with this patch, and no optimizations
currently take advantage of the new bool returned by do_add, do_multiply
and do_divide.
The following patch has been tested on i686-pc-linux-gnu with a full
"make bootstrap", all languages except Ada and treelang, and regression
tested with a top-level "make -k check" with no new failures.
Ok for mainline?
2003-04-20 Roger Sayle <roger at eyesopen dot com>
* real.c (do_add): Change to return a bool indicating that the
result of the operation may be inexact due to loss of precision.
(do_multiply): Likewise.
(do_divide): Likewise.
Index: real.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/real.c,v
retrieving revision 1.116
diff -c -3 -p -r1.116 real.c
*** real.c 2 Apr 2003 09:13:33 -0000 1.116
--- real.c 20 Apr 2003 16:24:37 -0000
*************** static bool div_significands PARAMS ((RE
*** 115,126 ****
const REAL_VALUE_TYPE *));
static void normalize PARAMS ((REAL_VALUE_TYPE *));
! static void do_add PARAMS ((REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
const REAL_VALUE_TYPE *, int));
! static void do_multiply PARAMS ((REAL_VALUE_TYPE *,
const REAL_VALUE_TYPE *,
const REAL_VALUE_TYPE *));
! static void do_divide PARAMS ((REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
const REAL_VALUE_TYPE *));
static int do_compare PARAMS ((const REAL_VALUE_TYPE *,
const REAL_VALUE_TYPE *, int));
--- 115,126 ----
const REAL_VALUE_TYPE *));
static void normalize PARAMS ((REAL_VALUE_TYPE *));
! static bool do_add PARAMS ((REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
const REAL_VALUE_TYPE *, int));
! static bool do_multiply PARAMS ((REAL_VALUE_TYPE *,
const REAL_VALUE_TYPE *,
const REAL_VALUE_TYPE *));
! static bool do_divide PARAMS ((REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
const REAL_VALUE_TYPE *));
static int do_compare PARAMS ((const REAL_VALUE_TYPE *,
const REAL_VALUE_TYPE *, int));
*************** normalize (r)
*** 563,571 ****
}
}
! /* Return R = A + (SUBTRACT_P ? -B : B). */
! static void
do_add (r, a, b, subtract_p)
REAL_VALUE_TYPE *r;
const REAL_VALUE_TYPE *a, *b;
--- 563,572 ----
}
}
! /* Calculate R = A + (SUBTRACT_P ? -B : B). Return true if the
! result may be inexact due to a loss of precision. */
! static bool
do_add (r, a, b, subtract_p)
REAL_VALUE_TYPE *r;
const REAL_VALUE_TYPE *a, *b;
*************** do_add (r, a, b, subtract_p)
*** 584,590 ****
case CLASS2 (rvc_zero, rvc_zero):
/* -0 + -0 = -0, -0 - +0 = -0; all other cases yield +0. */
get_zero (r, sign & !subtract_p);
! return;
case CLASS2 (rvc_zero, rvc_normal):
case CLASS2 (rvc_zero, rvc_inf):
--- 585,591 ----
case CLASS2 (rvc_zero, rvc_zero):
/* -0 + -0 = -0, -0 - +0 = -0; all other cases yield +0. */
get_zero (r, sign & !subtract_p);
! return false;
case CLASS2 (rvc_zero, rvc_normal):
case CLASS2 (rvc_zero, rvc_inf):
*************** do_add (r, a, b, subtract_p)
*** 598,604 ****
/* R + Inf = Inf. */
*r = *b;
r->sign = sign ^ subtract_p;
! return;
case CLASS2 (rvc_normal, rvc_zero):
case CLASS2 (rvc_inf, rvc_zero):
--- 599,605 ----
/* R + Inf = Inf. */
*r = *b;
r->sign = sign ^ subtract_p;
! return false;
case CLASS2 (rvc_normal, rvc_zero):
case CLASS2 (rvc_inf, rvc_zero):
*************** do_add (r, a, b, subtract_p)
*** 610,616 ****
case CLASS2 (rvc_inf, rvc_normal):
/* Inf + R = Inf. */
*r = *a;
! return;
case CLASS2 (rvc_inf, rvc_inf):
if (subtract_p)
--- 611,617 ----
case CLASS2 (rvc_inf, rvc_normal):
/* Inf + R = Inf. */
*r = *a;
! return false;
case CLASS2 (rvc_inf, rvc_inf):
if (subtract_p)
*************** do_add (r, a, b, subtract_p)
*** 619,625 ****
else
/* Inf + Inf = Inf. */
*r = *a;
! return;
case CLASS2 (rvc_normal, rvc_normal):
break;
--- 620,626 ----
else
/* Inf + Inf = Inf. */
*r = *a;
! return false;
case CLASS2 (rvc_normal, rvc_normal):
break;
*************** do_add (r, a, b, subtract_p)
*** 649,655 ****
{
*r = *a;
r->sign = sign;
! return;
}
inexact |= sticky_rshift_significand (&t, b, dexp);
--- 650,656 ----
{
*r = *a;
r->sign = sign;
! return true;
}
inexact |= sticky_rshift_significand (&t, b, dexp);
*************** do_add (r, a, b, subtract_p)
*** 680,686 ****
if (++exp > MAX_EXP)
{
get_inf (r, sign);
! return;
}
}
}
--- 681,687 ----
if (++exp > MAX_EXP)
{
get_inf (r, sign);
! return true;
}
}
}
*************** do_add (r, a, b, subtract_p)
*** 698,708 ****
r->sign = 0;
else
r->sig[0] |= inexact;
}
! /* Return R = A * B. */
! static void
do_multiply (r, a, b)
REAL_VALUE_TYPE *r;
const REAL_VALUE_TYPE *a, *b;
--- 699,711 ----
r->sign = 0;
else
r->sig[0] |= inexact;
+
+ return inexact;
}
! /* Calculate R = A * B. Return true if the result may be inexact. */
! static bool
do_multiply (r, a, b)
REAL_VALUE_TYPE *r;
const REAL_VALUE_TYPE *a, *b;
*************** do_multiply (r, a, b)
*** 710,715 ****
--- 713,719 ----
REAL_VALUE_TYPE u, t, *rr;
unsigned int i, j, k;
int sign = a->sign ^ b->sign;
+ bool inexact = false;
switch (CLASS2 (a->class, b->class))
{
*************** do_multiply (r, a, b)
*** 718,724 ****
case CLASS2 (rvc_normal, rvc_zero):
/* +-0 * ANY = 0 with appropriate sign. */
get_zero (r, sign);
! return;
case CLASS2 (rvc_zero, rvc_nan):
case CLASS2 (rvc_normal, rvc_nan):
--- 722,728 ----
case CLASS2 (rvc_normal, rvc_zero):
/* +-0 * ANY = 0 with appropriate sign. */
get_zero (r, sign);
! return false;
case CLASS2 (rvc_zero, rvc_nan):
case CLASS2 (rvc_normal, rvc_nan):
*************** do_multiply (r, a, b)
*** 727,733 ****
/* ANY * NaN = NaN. */
*r = *b;
r->sign = sign;
! return;
case CLASS2 (rvc_nan, rvc_zero):
case CLASS2 (rvc_nan, rvc_normal):
--- 731,737 ----
/* ANY * NaN = NaN. */
*r = *b;
r->sign = sign;
! return false;
case CLASS2 (rvc_nan, rvc_zero):
case CLASS2 (rvc_nan, rvc_normal):
*************** do_multiply (r, a, b)
*** 735,755 ****
/* NaN * ANY = NaN. */
*r = *a;
r->sign = sign;
! return;
case CLASS2 (rvc_zero, rvc_inf):
case CLASS2 (rvc_inf, rvc_zero):
/* 0 * Inf = NaN */
get_canonical_qnan (r, sign);
! return;
case CLASS2 (rvc_inf, rvc_inf):
case CLASS2 (rvc_normal, rvc_inf):
case CLASS2 (rvc_inf, rvc_normal):
/* Inf * Inf = Inf, R * Inf = Inf */
- overflow:
get_inf (r, sign);
! return;
case CLASS2 (rvc_normal, rvc_normal):
break;
--- 739,758 ----
/* NaN * ANY = NaN. */
*r = *a;
r->sign = sign;
! return false;
case CLASS2 (rvc_zero, rvc_inf):
case CLASS2 (rvc_inf, rvc_zero):
/* 0 * Inf = NaN */
get_canonical_qnan (r, sign);
! return false;
case CLASS2 (rvc_inf, rvc_inf):
case CLASS2 (rvc_normal, rvc_inf):
case CLASS2 (rvc_inf, rvc_normal):
/* Inf * Inf = Inf, R * Inf = Inf */
get_inf (r, sign);
! return false;
case CLASS2 (rvc_normal, rvc_normal):
break;
*************** do_multiply (r, a, b)
*** 799,808 ****
+ (b->exp - (1-j)*(HOST_BITS_PER_LONG/2)));
if (exp > MAX_EXP)
! goto overflow;
if (exp < -MAX_EXP)
! /* Would underflow to zero, which we shouldn't bother adding. */
! continue;
u.class = rvc_normal;
u.sign = 0;
--- 802,817 ----
+ (b->exp - (1-j)*(HOST_BITS_PER_LONG/2)));
if (exp > MAX_EXP)
! {
! get_inf (r, sign);
! return true;
! }
if (exp < -MAX_EXP)
! {
! /* Would underflow to zero, which we shouldn't bother adding. */
! inexact = true;
! continue;
! }
u.class = rvc_normal;
u.sign = 0;
*************** do_multiply (r, a, b)
*** 820,837 ****
}
normalize (&u);
! do_add (rr, rr, &u, 0);
}
}
rr->sign = sign;
if (rr != r)
*r = t;
}
! /* Return R = A / B. */
! static void
do_divide (r, a, b)
REAL_VALUE_TYPE *r;
const REAL_VALUE_TYPE *a, *b;
--- 829,848 ----
}
normalize (&u);
! inexact |= do_add (rr, rr, &u, 0);
}
}
rr->sign = sign;
if (rr != r)
*r = t;
+
+ return inexact;
}
! /* Calculate R = A / B. Return true if the result may be inexact. */
! static bool
do_divide (r, a, b)
REAL_VALUE_TYPE *r;
const REAL_VALUE_TYPE *a, *b;
*************** do_divide (r, a, b)
*** 847,869 ****
case CLASS2 (rvc_inf, rvc_inf):
/* Inf / Inf = NaN. */
get_canonical_qnan (r, sign);
! return;
case CLASS2 (rvc_zero, rvc_normal):
case CLASS2 (rvc_zero, rvc_inf):
/* 0 / ANY = 0. */
case CLASS2 (rvc_normal, rvc_inf):
/* R / Inf = 0. */
- underflow:
get_zero (r, sign);
! return;
case CLASS2 (rvc_normal, rvc_zero):
/* R / 0 = Inf. */
case CLASS2 (rvc_inf, rvc_zero):
/* Inf / 0 = Inf. */
get_inf (r, sign);
! return;
case CLASS2 (rvc_zero, rvc_nan):
case CLASS2 (rvc_normal, rvc_nan):
--- 858,879 ----
case CLASS2 (rvc_inf, rvc_inf):
/* Inf / Inf = NaN. */
get_canonical_qnan (r, sign);
! return false;
case CLASS2 (rvc_zero, rvc_normal):
case CLASS2 (rvc_zero, rvc_inf):
/* 0 / ANY = 0. */
case CLASS2 (rvc_normal, rvc_inf):
/* R / Inf = 0. */
get_zero (r, sign);
! return false;
case CLASS2 (rvc_normal, rvc_zero):
/* R / 0 = Inf. */
case CLASS2 (rvc_inf, rvc_zero):
/* Inf / 0 = Inf. */
get_inf (r, sign);
! return false;
case CLASS2 (rvc_zero, rvc_nan):
case CLASS2 (rvc_normal, rvc_nan):
*************** do_divide (r, a, b)
*** 872,878 ****
/* ANY / NaN = NaN. */
*r = *b;
r->sign = sign;
! return;
case CLASS2 (rvc_nan, rvc_zero):
case CLASS2 (rvc_nan, rvc_normal):
--- 882,888 ----
/* ANY / NaN = NaN. */
*r = *b;
r->sign = sign;
! return false;
case CLASS2 (rvc_nan, rvc_zero):
case CLASS2 (rvc_nan, rvc_normal):
*************** do_divide (r, a, b)
*** 880,892 ****
/* NaN / ANY = NaN. */
*r = *a;
r->sign = sign;
! return;
case CLASS2 (rvc_inf, rvc_normal):
/* Inf / R = Inf. */
- overflow:
get_inf (r, sign);
! return;
case CLASS2 (rvc_normal, rvc_normal):
break;
--- 890,901 ----
/* NaN / ANY = NaN. */
*r = *a;
r->sign = sign;
! return false;
case CLASS2 (rvc_inf, rvc_normal):
/* Inf / R = Inf. */
get_inf (r, sign);
! return false;
case CLASS2 (rvc_normal, rvc_normal):
break;
*************** do_divide (r, a, b)
*** 905,913 ****
exp = a->exp - b->exp + 1;
if (exp > MAX_EXP)
! goto overflow;
if (exp < -MAX_EXP)
! goto underflow;
rr->exp = exp;
inexact = div_significands (rr, a, b);
--- 914,928 ----
exp = a->exp - b->exp + 1;
if (exp > MAX_EXP)
! {
! get_inf (r, sign);
! return true;
! }
if (exp < -MAX_EXP)
! {
! get_zero (r, sign);
! return true;
! }
rr->exp = exp;
inexact = div_significands (rr, a, b);
*************** do_divide (r, a, b)
*** 918,923 ****
--- 933,940 ----
if (rr != r)
*r = t;
+
+ return inexact;
}
/* Return a tri-state comparison of A vs B. Return NAN_RESULT if
Roger
--
Roger Sayle, E-mail: roger at eyesopen dot com
OpenEye Scientific Software, WWW: http://www.eyesopen.com/
Suite 1107, 3600 Cerrillos Road, Tel: (+1) 505-473-7385
Santa Fe, New Mexico, 87507. Fax: (+1) 505-473-0833