[glibc] powerpc: Update acosf ulps

[Adhemerval Zanella Netto]

On 02/01/25 10:32, Paul Zimmermann wrote:
>        Hi Alexander,
> 
>>>> Isn't the new acosf implementation supposed to be correctly rounded?
>>>
>>> yes, but due to https://gcc.gnu.org/bugzilla/show_bug.cgi?id=57245,
>>> the new implementation might be miscompiled by gcc < 12.
>>
>> can you kindly point me where the new acosf employs problematic conversions?
>> I see a couple of instances where the 'f' suffix is missing, producing a
>> double literal where a float could work as well:
>>
>>  * in as_special there is
>>
>>   const float pih = 0x1.921fb6p+1;
> 
> should be changed to 0x1.921fb6p+1f for more clarity
> 
>>  * and in __ieee754_acosf, in 0x1p-25:
>>
>>       if (t == 0x328885a3u)
>>        return 0x1.921fb6p+0f + 0x1p-25;
>>       if (t == 0x39826222u)
>>        return 0x1.920f6ap+0f + 0x1p-25;
> 
> both 0x1p-25 should be 0x1p-25f for more clarity
> 
>> but I don't see any instance where a double-precision constant is converted
>> to float such that result depends on rounding mode.
> 
> indeed, this does not change the issue with gcc <= 11
> 
>> I'm guessing the '+ 0x1p-25' trick above is a workaround for this exact bug?
> 
> no, this is to get correct rounding depending on the rounding mode.
> 

This seems to show no regression on gcc 11 and gcc 13 on x86_64, i686, aarch64,
and powerpc64le:

diff --git a/sysdeps/ieee754/flt-32/e_acosf.c b/sysdeps/ieee754/flt-32/e_acosf.c
index cba01221dc..35e241f6c4 100644
--- a/sysdeps/ieee754/flt-32/e_acosf.c
+++ b/sysdeps/ieee754/flt-32/e_acosf.c
@@ -28,6 +28,7 @@ SOFTWARE.
 #include <math.h>
 #include <math_private.h>
 #include <libm-alias-finite.h>
+#include <math-barriers.h>
 #include "math_config.h"

 static __attribute__ ((noinline)) float
@@ -66,7 +67,7 @@ poly12 (double z, const double *c)
 float
 __ieee754_acosf (float x)
 {
-  const double pi2 = 0x1.921fb54442d18p+0;
+  double pi2 = 0x1.921fb54442d18p+0;
   static const double o[] = { 0, 0x1.921fb54442d18p+1 };
   double xs = x;
   double r;
@@ -87,7 +88,10 @@ __ieee754_acosf (float x)
        };
       /* Avoid spurious underflow exception.  */
       if (__glibc_unlikely (ax <= 0x40000000u)) /* |x| < 2^-63 */
-       return (float) pi2;
+       /* GCC <= 11 wrongly assumes the rounding is to nearest and
+          performs a constant folding here:
+          https://gcc.gnu.org/bugzilla/show_bug.cgi?id=57245 */
+       return math_opt_barrier (pi2);
       double z = xs;
       double z2 = z * z;
       double z4 = z2 * z2;