From 95c26233a1c7dd73a852a4f3bffd2d2eeba5e2fd Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Fri, 13 Feb 2015 16:20:36 +0000 Subject: [PATCH] Fix powerpc software sqrtf (bug 17967). Similarly to sqrt in , the powerpc sqrtf implementation for when _ARCH_PPCSQ is not defined also relies on a * b + c being contracted into a fused multiply-add. Although this contraction is not explicitly disabled for e_sqrtf.c, it still seems appropriate to make the file explicit about its requirements by using __builtin_fmaf; this patch does so. Furthermore, it turns out that doing so fixes the observed inaccuracy and missing exceptions (that is, that without explicit __builtin_fmaf usage, it was not being compiled as intended). Tested for powerpc32 (hard float). [BZ #17967] * sysdeps/powerpc/fpu/e_sqrtf.c (__slow_ieee754_sqrtf): Use __builtin_fmaf instead of relying on contraction of a * b + c. --- ChangeLog | 6 ++++++ NEWS | 2 +- sysdeps/powerpc/fpu/e_sqrtf.c | 24 +++++++++++++----------- 3 files changed, 20 insertions(+), 12 deletions(-) diff --git a/ChangeLog b/ChangeLog index 821da1dd88..50adb3cbb3 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,9 @@ +2015-02-13 Joseph Myers + + [BZ #17967] + * sysdeps/powerpc/fpu/e_sqrtf.c (__slow_ieee754_sqrtf): Use + __builtin_fmaf instead of relying on contraction of a * b + c. + 2015-02-12 J William Piggott [BZ #17969] diff --git a/NEWS b/NEWS index f6f34e8811..b95f5e6d25 100644 --- a/NEWS +++ b/NEWS @@ -10,7 +10,7 @@ Version 2.22 * The following bugs are resolved with this release: 4719, 15467, 15790, 16560, 17912, 17932, 17944, 17949, 17964, 17965, - 17969. + 17967, 17969. Version 2.21 diff --git a/sysdeps/powerpc/fpu/e_sqrtf.c b/sysdeps/powerpc/fpu/e_sqrtf.c index 034b6f58ca..a684cf977a 100644 --- a/sysdeps/powerpc/fpu/e_sqrtf.c +++ b/sysdeps/powerpc/fpu/e_sqrtf.c @@ -87,26 +87,28 @@ __slow_ieee754_sqrtf (float x) /* Here we have three Newton-Raphson iterations each of a division and a square root and the remainder of the argument reduction, all interleaved. */ - sd = -(sg * sg - sx); + sd = -__builtin_fmaf (sg, sg, -sx); fsgi = (xi + 0x40000000) >> 1 & 0x7f800000; sy2 = sy + sy; - sg = sy * sd + sg; /* 16-bit approximation to sqrt(sx). */ - e = -(sy * sg - almost_half); + sg = __builtin_fmaf (sy, sd, sg); /* 16-bit approximation to + sqrt(sx). */ + e = -__builtin_fmaf (sy, sg, -almost_half); SET_FLOAT_WORD (fsg, fsgi); - sd = -(sg * sg - sx); - sy = sy + e * sy2; + sd = -__builtin_fmaf (sg, sg, -sx); + sy = __builtin_fmaf (e, sy2, sy); if ((xi & 0x7f800000) == 0) goto denorm; shx = sx * fsg; - sg = sg + sy * sd; /* 32-bit approximation to sqrt(sx), - but perhaps rounded incorrectly. */ + sg = __builtin_fmaf (sy, sd, sg); /* 32-bit approximation to + sqrt(sx), but perhaps + rounded incorrectly. */ sy2 = sy + sy; g = sg * fsg; - e = -(sy * sg - almost_half); - d = -(g * sg - shx); - sy = sy + e * sy2; + e = -__builtin_fmaf (sy, sg, -almost_half); + d = -__builtin_fmaf (g, sg, -shx); + sy = __builtin_fmaf (e, sy2, sy); fesetenv_register (fe); - return g + sy * d; + return __builtin_fmaf (sy, d, g); denorm: /* For denormalised numbers, we normalise, calculate the square root, and return an adjusted result. */ -- 2.11.4.GIT