[PATCH 1/3] Optimize trunc() and truncf().

[PATCH 1/3] Optimize trunc() and truncf().

[Matt Turner]

By creating a mask of non-fractional bits from the exponent.
---

Joseph suggested an SSE 4.1 implementation of trunc/truncf, so I thought
now would be a good time to send these patches even sans benchmarking data.

I do not believe the other two generic trunc* implementations (ldbl-128,
dbl-64) would benefit from the change made in this patch

Suggestions for ChangeLog entries welcome. Guidance for generating benchmark
data requested.

 sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c | 32 +++++++++++-----------------
 sysdeps/ieee754/flt-32/s_truncf.c            | 32 +++++++++++-----------------
 2 files changed, 26 insertions(+), 38 deletions(-)

diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c
index 81ac55e..e4cba3b 100644
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c
+++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c
@@ -1,7 +1,6 @@
 /* Truncate argument to nearest integral value not larger than the argument.
    Copyright (C) 1997-2016 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
-   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
    The GNU C Library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
@@ -21,31 +20,26 @@
 
 #include <math_private.h>
 
+static int64_t
+max (int64_t x, int64_t y)
+{
+  return x > y ? x : y;
+}
 
 double
 __trunc (double x)
 {
-  int64_t i0, j0;
-  int64_t sx;
+  int64_t i0;
 
   EXTRACT_WORDS64 (i0, x);
-  sx = i0 & UINT64_C(0x8000000000000000);
-  j0 = ((i0 >> 52) & 0x7ff) - 0x3ff;
-  if (j0 < 52)
-    {
-      if (j0 < 0)
-	/* The magnitude of the number is < 1 so the result is +-0.  */
-	INSERT_WORDS64 (x, sx);
-      else
-	INSERT_WORDS64 (x, sx | (i0 & ~(UINT64_C(0x000fffffffffffff) >> j0)));
-    }
-  else
-    {
-      if (j0 == 0x400)
-	/* x is inf or NaN.  */
-	return x + x;
-    }
+  int64_t exp = (i0 >> 52) & 0x7ff;
+  int64_t mask = UINT64_C(-1) << max(52 - (exp - 1023), 0);
+
+  if (exp < 1023)
+    mask = UINT64_C(0x8000000000000000);
 
+  i0 &= mask;
+  INSERT_WORDS64(x, i0);
   return x;
 }
 weak_alias (__trunc, trunc)
diff --git a/sysdeps/ieee754/flt-32/s_truncf.c b/sysdeps/ieee754/flt-32/s_truncf.c
index 43d35c7..67fdcc8 100644
--- a/sysdeps/ieee754/flt-32/s_truncf.c
+++ b/sysdeps/ieee754/flt-32/s_truncf.c
@@ -1,7 +1,6 @@
 /* Truncate argument to nearest integral value not larger than the argument.
    Copyright (C) 1997-2016 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
-   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
    The GNU C Library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
@@ -21,31 +20,26 @@
 
 #include <math_private.h>
 
+static int32_t
+max (int32_t x, int32_t y)
+{
+  return x > y ? x : y;
+}
 
 float
 __truncf (float x)
 {
-  int32_t i0, j0;
-  int sx;
+  int32_t i0;
 
   GET_FLOAT_WORD (i0, x);
-  sx = i0 & 0x80000000;
-  j0 = ((i0 >> 23) & 0xff) - 0x7f;
-  if (j0 < 23)
-    {
-      if (j0 < 0)
-	/* The magnitude of the number is < 1 so the result is +-0.  */
-	SET_FLOAT_WORD (x, sx);
-      else
-	SET_FLOAT_WORD (x, sx | (i0 & ~(0x007fffff >> j0)));
-    }
-  else
-    {
-      if (j0 == 0x80)
-	/* x is inf or NaN.  */
-	return x + x;
-    }
+  int32_t exp = (i0 >> 23) & 0xff;
+  int32_t mask = ~0u << max(23 - (exp - 127), 0);
+
+  if (exp < 127)
+    mask = 0x80000000;
 
+  i0 &= mask;
+  SET_FLOAT_WORD (x, i0);
   return x;
 }
 weak_alias (__truncf, truncf)
-- 
2.7.3

[PATCH 2/3] Add more tests for trunc.

[Matt Turner]

Some implementations of truncf(x) add 2**23 to x with the rounding mode
set to truncate before subtracting 2**23 to produce the truncated
result. Unfortunately this doesn't work for odd values of x greater than
2**23, since the result of the addition is inexact and the parity bit is
lost by the rounding step. The same problems apply to implementations of
trunc on other floating-point data types.

2016-05-24  Matt Turner  <mattst88@gmail.com>

        * math/libm-test.inc (trunc_test_data): Add tests for odd
        values.
---
 math/libm-test.inc | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/math/libm-test.inc b/math/libm-test.inc
index f1ba7dd..84c119a 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -11413,18 +11413,23 @@ static const struct test_f_f_data trunc_test_data[] =
     TEST_f_f (trunc, 1, 1, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, 2, 2, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, 0x1p23, 0x1p23, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
+    TEST_f_f (trunc, 0x1p23 + 1, 0x1p23 + 1, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, 0x1p24, 0x1p24, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, 0x1p25, 0x1p25, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, 0x1p52, 0x1p52, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
+    TEST_f_f (trunc, 0x1p52 + 1, 0x1p52 + 1, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, 0x1p53, 0x1p53, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, 0x1p54, 0x1p54, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, 0x1p63, 0x1p63, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
+    TEST_f_f (trunc, 0x1p63 + 1, 0x1p63 + 1, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, 0x1p64, 0x1p64, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, 0x1p65, 0x1p65, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, 0x1p105, 0x1p105, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
+    TEST_f_f (trunc, 0x1p105 + 1, 0x1p105 + 1, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, 0x1p106, 0x1p106, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, 0x1p107, 0x1p107, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, 0x1p112, 0x1p112, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
+    TEST_f_f (trunc, 0x1p112 + 1, 0x1p112 + 1, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, 0x1p113, 0x1p113, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, 0x1p114, 0x1p114, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, max_value, max_value, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
@@ -11433,18 +11438,23 @@ static const struct test_f_f_data trunc_test_data[] =
     TEST_f_f (trunc, -1, -1, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, -2, -2, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, -0x1p23, -0x1p23, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
+    TEST_f_f (trunc, -0x1p23 - 1, -0x1p23 - 1, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, -0x1p24, -0x1p24, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, -0x1p25, -0x1p25, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, -0x1p52, -0x1p52, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
+    TEST_f_f (trunc, -0x1p52 - 1, -0x1p52 - 1, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, -0x1p53, -0x1p53, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, -0x1p54, -0x1p54, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, -0x1p63, -0x1p63, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
+    TEST_f_f (trunc, -0x1p63 - 1, -0x1p63 - 1, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, -0x1p64, -0x1p64, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, -0x1p65, -0x1p65, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, -0x1p105, -0x1p105, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
+    TEST_f_f (trunc, -0x1p105 - 1, -0x1p105 - 1, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, -0x1p106, -0x1p106, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, -0x1p107, -0x1p107, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, -0x1p112, -0x1p112, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
+    TEST_f_f (trunc, -0x1p112 - 1, -0x1p112 - 1, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, -0x1p113, -0x1p113, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, -0x1p114, -0x1p114, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (trunc, -max_value, -max_value, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
-- 
2.7.3

[PATCH 3/3] alpha: Remove trunc() and truncf() implementations.

[Matt Turner]

They returned incorrect results for odd values greater than 2**52 and
2**23 respectively; and also raised incorrect inexact exceptions.

2016-05-24  Matt Turner  <mattst88@gmail.com>

    * sysdeps/alpha/fpu/s_trunc.c: Remove file.
    * sysdeps/alpha/fpu/s_truncf.c: Likewise.
---
 sysdeps/alpha/fpu/s_trunc.c  | 52 --------------------------------------------
 sysdeps/alpha/fpu/s_truncf.c | 44 -------------------------------------
 2 files changed, 96 deletions(-)
 delete mode 100644 sysdeps/alpha/fpu/s_trunc.c
 delete mode 100644 sysdeps/alpha/fpu/s_truncf.c

diff --git a/sysdeps/alpha/fpu/s_trunc.c b/sysdeps/alpha/fpu/s_trunc.c
deleted file mode 100644
index 16cb114..0000000
--- a/sysdeps/alpha/fpu/s_trunc.c
+++ /dev/null
@@ -1,52 +0,0 @@
-/* Copyright (C) 2007-2016 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-   Contributed by Richard Henderson.
-
-   The GNU C Library is free software; you can redistribute it and/or
-   modify it under the terms of the GNU Lesser General Public
-   License as published by the Free Software Foundation; either
-   version 2.1 of the License, or (at your option) any later version.
-
-   The GNU C Library 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
-   Lesser General Public License for more details.
-
-   You should have received a copy of the GNU Lesser General Public
-   License along with the GNU C Library.  If not, see
-   <http://www.gnu.org/licenses/>.  */
-
-#include <math.h>
-#include <math_ldbl_opt.h>
-
-
-/* Use the chopped rounding mode conversion instructions to implement trunc. */
-
-double
-__trunc (double x)
-{
-  double two52 = copysign (0x1.0p52, x);
-  double r, tmp;
-
-  __asm (
-#ifdef _IEEE_FP_INEXACT
-	 "addt/suic %2, %3, %1\n\tsubt/suic %1, %3, %0"
-#else
-	 "addt/suc %2, %3, %1\n\tsubt/suc %1, %3, %0"
-#endif
-	 : "=&f"(r), "=&f"(tmp)
-	 : "f"(x), "f"(two52));
-
-  /* trunc(-0) == -0, and in general we'll always have the same
-     sign as our input.  */
-  return copysign (r, x);
-}
-
-weak_alias (__trunc, trunc)
-#ifdef NO_LONG_DOUBLE
-strong_alias (__trunc, __truncl)
-weak_alias (__trunc, truncl)
-#endif
-#if LONG_DOUBLE_COMPAT(libm, GLIBC_2_1)
-compat_symbol (libm, __trunc, truncl, GLIBC_2_1);
-#endif
diff --git a/sysdeps/alpha/fpu/s_truncf.c b/sysdeps/alpha/fpu/s_truncf.c
deleted file mode 100644
index 2290f28..0000000
--- a/sysdeps/alpha/fpu/s_truncf.c
+++ /dev/null
@@ -1,44 +0,0 @@
-/* Copyright (C) 2007-2016 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-   Contributed by Richard Henderson.
-
-   The GNU C Library is free software; you can redistribute it and/or
-   modify it under the terms of the GNU Lesser General Public
-   License as published by the Free Software Foundation; either
-   version 2.1 of the License, or (at your option) any later version.
-
-   The GNU C Library 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
-   Lesser General Public License for more details.
-
-   You should have received a copy of the GNU Lesser General Public
-   License along with the GNU C Library.  If not, see
-   <http://www.gnu.org/licenses/>.  */
-
-#include <math.h>
-
-
-/* Use the chopped rounding mode conversion instructions to implement trunc. */
-
-float
-__truncf (float x)
-{
-  float two23 = copysignf (0x1.0p23, x);
-  float r, tmp;
-
-  __asm (
-#ifdef _IEEE_FP_INEXACT
-	 "adds/suic %2, %3, %1\n\tsubs/suic %1, %3, %0"
-#else
-	 "adds/suc %2, %3, %1\n\tsubs/suc %1, %3, %0"
-#endif
-	 : "=&f"(r), "=&f"(tmp)
-	 : "f"(x), "f"(two23));
-
-  /* trunc(-0) == -0, and in general we'll always have the same
-     sign as our input.  */
-  return copysignf (r, x);
-}
-
-weak_alias (__truncf, truncf)
-- 
2.7.3

Re: [PATCH 1/3] Optimize trunc() and truncf().

[Florian Weimer]

On 05/25/2016 03:32 AM, Matt Turner wrote:
 > +  int64_t exp = (i0 >> 52) & 0x7ff;
 > +  int64_t mask = UINT64_C(-1) << max(52 - (exp - 1023), 0);

I think it will help GCC if you make exp an int and change the max 
function accordingly.

You should use spaces in front of the parenthesis in a function call, 
and also for most macro calls.

Thanks,
Florian

Re: [PATCH 1/3] Optimize trunc() and truncf().

[Joseph Myers]

On Tue, 24 May 2016, Matt Turner wrote:

> +  int64_t exp = (i0 >> 52) & 0x7ff;
> +  int64_t mask = UINT64_C(-1) << max(52 - (exp - 1023), 0);

This can involve undefined behavior (shift by more than 63), so ...

> +  if (exp < 1023)
> +    mask = UINT64_C(0x8000000000000000);

... the shift should only be done at all if exp is large enough.

Likewise for truncf.

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: [PATCH 2/3] Add more tests for trunc.

[Joseph Myers]

On Tue, 24 May 2016, Matt Turner wrote:

> Some implementations of truncf(x) add 2**23 to x with the rounding mode
> set to truncate before subtracting 2**23 to produce the truncated
> result. Unfortunately this doesn't work for odd values of x greater than
> 2**23, since the result of the addition is inexact and the parity bit is
> lost by the rounding step. The same problems apply to implementations of
> trunc on other floating-point data types.

I think the same tests should be added for all of ceil, floor, round, 
trunc if not already present.

> +    TEST_f_f (trunc, 0x1p63 + 1, 0x1p63 + 1, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
>      TEST_f_f (trunc, 0x1p64, 0x1p64, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
>      TEST_f_f (trunc, 0x1p65, 0x1p65, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
>      TEST_f_f (trunc, 0x1p105, 0x1p105, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
> +    TEST_f_f (trunc, 0x1p105 + 1, 0x1p105 + 1, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
>      TEST_f_f (trunc, 0x1p106, 0x1p106, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
>      TEST_f_f (trunc, 0x1p107, 0x1p107, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
>      TEST_f_f (trunc, 0x1p112, 0x1p112, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
> +    TEST_f_f (trunc, 0x1p112 + 1, 0x1p112 + 1, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),

These need to suffix the constants with 'L' so that they actually test 
what's intended for long double.  And doing arithmetic like that in 
libm-test.inc makes manipulation to change suffixes, as in 
<https://sourceware.org/ml/libc-alpha/2016-05/msg00474.html>, harder.  So 
I think it's better to write all the constants explicitly without 
arithmetic (0x1.000002p23L, 0x1.0000000000001p52L, 
0x1.0000000000000002p63L, 0x1.000000000000000000000000008p105L, 
0x1.0000000000000000000000000001p112L).

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: [PATCH 3/3] alpha: Remove trunc() and truncf() implementations.

[Joseph Myers]

On Tue, 24 May 2016, Matt Turner wrote:

> They returned incorrect results for odd values greater than 2**52 and
> 2**23 respectively; and also raised incorrect inexact exceptions.

Incorrect results means a bug that's not currently filed in Bugzilla, and 
as a bug user-visible in releases it should be filed there before fixing 
(and then [BZ #N] included in the ChangeLog entry and the bug resolved as 
FIXED with milestone set accordingly when the fix goes in, so that the bug 
appears in the list of fixed bugs generated for NEWS).

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: [PATCH 1/3] Optimize trunc() and truncf().

[Joseph Myers]

On Tue, 24 May 2016, Matt Turner wrote:

> -  else
> -    {
> -      if (j0 == 0x400)
> -	/* x is inf or NaN.  */
> -	return x + x;
> -    }

You're removing too much code here.  You need to keep this addition in the 
NaN case to ensure that signaling NaNs raise "invalid" and get converted 
to quiet NaNs in the process.

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: [PATCH 1/3] Optimize trunc() and truncf().

[Mike Frysinger]

[-- Attachment #1: Type: text/plain, Size: 462 bytes --]

On 25 May 2016 15:36, Joseph Myers wrote:
> On Tue, 24 May 2016, Matt Turner wrote:
> > -  else
> > -    {
> > -      if (j0 == 0x400)
> > -	/* x is inf or NaN.  */
> > -	return x + x;
> > -    }
> 
> You're removing too much code here.  You need to keep this addition in the 
> NaN case to ensure that signaling NaNs raise "invalid" and get converted 
> to quiet NaNs in the process.

does this mean current test coverage is insufficient ?
-mike

[-- Attachment #2: Digital signature --]
[-- Type: application/pgp-signature, Size: 819 bytes --]

Re: [PATCH 1/3] Optimize trunc() and truncf().

[Joseph Myers]

On Wed, 25 May 2016, Mike Frysinger wrote:

> On 25 May 2016 15:36, Joseph Myers wrote:
> > On Tue, 24 May 2016, Matt Turner wrote:
> > > -  else
> > > -    {
> > > -      if (j0 == 0x400)
> > > -	/* x is inf or NaN.  */
> > > -	return x + x;
> > > -    }
> > 
> > You're removing too much code here.  You need to keep this addition in the 
> > NaN case to ensure that signaling NaNs raise "invalid" and get converted 
> > to quiet NaNs in the process.
> 
> does this mean current test coverage is insufficient ?

Yes.  I intend to add sNaN support to libm-test.inc at some point if noone 
else gets there first; no doubt this would then show up lots of bugs in 
existing code.  Cf. 
<https://sourceware.org/ml/libc-ports/2013-04/msg00008.html>, but that 
would need completely reworking now for current libm-test.

-- 
Joseph S. Myers
joseph@codesourcery.com