# HG changeset patch # User Bruno Haible # Date 1192233829 -7200 # Node ID 5437ef3873b58460387dc9f2c3f1c5d2f934dce4 # Parent cb8282d0ed8ceb77f44737877fe164d59b0ada3a Fix incorrect rounding of trunc, truncf, truncl in some cases. Add a new test. diff --git a/ChangeLog b/ChangeLog --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,10 @@ 2007-10-12 Bruno Haible + * lib/trunc.c (FUNC): Avoid rounding errors for values near a power + of 2. + * tests/test-truncf2.c: New file. + * modules/truncf-tests: (Files, Depends-on, Makefile.am): Add new test. + * tests/test-truncf1.c: Renamed from tests/test-truncf.c. * modules/truncf-tests: Update. diff --git a/lib/trunc.c b/lib/trunc.c --- a/lib/trunc.c +++ b/lib/trunc.c @@ -65,21 +65,29 @@ if (z > L_(0.0)) { - /* Round to the next integer (nearest or up or down, doesn't matter). */ - z += TWO_MANT_DIG; - z -= TWO_MANT_DIG; - /* Enforce rounding down. */ - if (z > y) - z -= L_(1.0); + /* Avoid rounding errors for values near 2^k, where k >= MANT_DIG-1. */ + if (z < TWO_MANT_DIG) + { + /* Round to the next integer (nearest or up or down, doesn't matter). */ + z += TWO_MANT_DIG; + z -= TWO_MANT_DIG; + /* Enforce rounding down. */ + if (z > y) + z -= L_(1.0); + } } else if (z < L_(0.0)) { - /* Round to the next integer (nearest or up or down, doesn't matter). */ - z -= TWO_MANT_DIG; - z += TWO_MANT_DIG; - /* Enforce rounding up. */ - if (z < y) - z += L_(1.0); + /* Avoid rounding errors for values near -2^k, where k >= MANT_DIG-1. */ + if (z > - TWO_MANT_DIG) + { + /* Round to the next integer (nearest or up or down, doesn't matter). */ + z -= TWO_MANT_DIG; + z += TWO_MANT_DIG; + /* Enforce rounding up. */ + if (z < y) + z += L_(1.0); + } } return z; } diff --git a/modules/truncf-tests b/modules/truncf-tests --- a/modules/truncf-tests +++ b/modules/truncf-tests @@ -1,14 +1,20 @@ Files: tests/test-truncf1.c +tests/test-truncf2.c Depends-on: +float +isnanf-nolibm +stdbool +stdint configure.ac: Makefile.am: -TESTS += test-truncf1 -check_PROGRAMS += test-truncf1 +TESTS += test-truncf1 test-truncf2 +check_PROGRAMS += test-truncf1 test-truncf2 test_truncf1_LDADD = $(LDADD) @TRUNCF_LIBM@ +test_truncf2_LDADD = $(LDADD) @TRUNCF_LIBM@ License: LGPL diff --git a/tests/test-truncf2.c b/tests/test-truncf2.c new file mode 100644 --- /dev/null +++ b/tests/test-truncf2.c @@ -0,0 +1,166 @@ +/* Test of rounding towards zero. + Copyright (C) 2007 Free Software Foundation, Inc. + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . */ + +/* Written by Bruno Haible , 2007. */ + +#include + +#include + +#include +#include +#include +#include +#include + +#include "isnanf.h" + +#define ASSERT(expr) \ + do \ + { \ + if (!(expr)) \ + { \ + fprintf (stderr, "%s:%d: assertion failed\n", __FILE__, __LINE__); \ + abort (); \ + } \ + } \ + while (0) + + +/* The reference implementation, taken from lib/trunc.c. */ + +#define DOUBLE float +#define MANT_DIG FLT_MANT_DIG +#define L_(literal) literal##f + +/* 2^(MANT_DIG-1). */ +static const DOUBLE TWO_MANT_DIG = + /* Assume MANT_DIG <= 5 * 31. + Use the identity + n = floor(n/5) + floor((n+1)/5) + ... + floor((n+4)/5). */ + (DOUBLE) (1U << ((MANT_DIG - 1) / 5)) + * (DOUBLE) (1U << ((MANT_DIG - 1 + 1) / 5)) + * (DOUBLE) (1U << ((MANT_DIG - 1 + 2) / 5)) + * (DOUBLE) (1U << ((MANT_DIG - 1 + 3) / 5)) + * (DOUBLE) (1U << ((MANT_DIG - 1 + 4) / 5)); + +DOUBLE +truncf_reference (DOUBLE x) +{ + /* The use of 'volatile' guarantees that excess precision bits are dropped + at each addition step and before the following comparison at the caller's + site. It is necessary on x86 systems where double-floats are not IEEE + compliant by default, to avoid that the results become platform and compiler + option dependent. 'volatile' is a portable alternative to gcc's + -ffloat-store option. */ + volatile DOUBLE y = x; + volatile DOUBLE z = y; + + if (z > L_(0.0)) + { + /* Avoid rounding errors for values near 2^k, where k >= MANT_DIG-1. */ + if (z < TWO_MANT_DIG) + { + /* Round to the next integer (nearest or up or down, doesn't matter). */ + z += TWO_MANT_DIG; + z -= TWO_MANT_DIG; + /* Enforce rounding down. */ + if (z > y) + z -= L_(1.0); + } + } + else if (z < L_(0.0)) + { + /* Avoid rounding errors for values near -2^k, where k >= MANT_DIG-1. */ + if (z > - TWO_MANT_DIG) + { + /* Round to the next integer (nearest or up or down, doesn't matter). */ + z -= TWO_MANT_DIG; + z += TWO_MANT_DIG; + /* Enforce rounding up. */ + if (z < y) + z += L_(1.0); + } + } + return z; +} + + +/* Test for equality. */ +static int +equal (DOUBLE x, DOUBLE y) +{ + return (isnanf (x) ? isnanf (y) : x == y); +} + +/* Test whether the result for a given argument is correct. */ +static bool +correct_result_p (DOUBLE x, DOUBLE result) +{ + return + (x >= 0 + ? (x < 1 ? result == L_(0.0) : + x - 1 < x ? result <= x && result > x - 1 : + equal (result, x)) + : (x > -1 ? result == L_(0.0) : + x + 1 > x ? result >= x && result < x + 1 : + equal (result, x))); +} + +/* Test the function for a given argument. */ +static int +check (float x) +{ + /* If the reference implementation is incorrect, bail out immediately. */ + float reference = truncf_reference (x); + ASSERT (correct_result_p (x, reference)); + /* If the actual implementation is wrong, return an error code. */ + { + float result = truncf (x); + if (correct_result_p (x, result)) + return 0; + else + { + fprintf (stderr, "truncf %g(%a) = %g(%a) or %g(%a)?\n", + x, x, reference, reference, result, result); + return 1; + } + } +} + +#define NUM_HIGHBITS 12 +#define NUM_LOWBITS 4 + +int +main () +{ + unsigned int highbits; + unsigned int lowbits; + int error = 0; + for (highbits = 0; highbits < (1 << NUM_HIGHBITS); highbits++) + for (lowbits = 0; lowbits < (1 << NUM_LOWBITS); lowbits++) + { + /* Combine highbits and lowbits into a floating-point number, + sign-extending the lowbits to 32-NUM_HIGHBITS bits. */ + union { float f; uint32_t i; } janus; + janus.i = ((uint32_t) highbits << (32 - NUM_HIGHBITS)) + | ((uint32_t) ((int32_t) ((uint32_t) lowbits << (32 - NUM_LOWBITS)) + >> (32 - NUM_LOWBITS - NUM_HIGHBITS)) + >> NUM_HIGHBITS); + error |= check (janus.f); + } + return (error ? 1 : 0); +}