Optimization bug with floating-point?

All,

There seems to an optimization bug with clang on

% uname -a
FreeBSD mobile 13.0-CURRENT FreeBSD 13.0-CURRENT r344653 MOBILE  i386

IOW, if you do numerica work on i386, you may want to check your
results.

The program demonstrating the issue is at the end of this email.

gcc8 --version
gcc8 (FreeBSD Ports Collection) 8.3.0

gcc8 -fno-builtin -o z a.c -lm && ./z
gcc8 -O -fno-builtin -o z a.c -lm && ./z
gcc8 -O2 -fno-builtin -o z a.c -lm && ./z
gcc8 -O3 -fno-builtin -o z a.c -lm && ./z

Max ULP: 2.297073
Count: 0           (# of ULP that exceed 21)


cc --version
FreeBSD clang version 7.0.1 (tags/RELEASE_701/final 349250)
(based on LLVM 7.0.1)
Target: i386-unknown-freebsd13.0

cc -fno-builtin -o z a.c -lm && ./z
Max ULP: 2.297073
Count: 0

cc -O -fno-builtin -o z a.c -lm && ./z
cc -O2 -fno-builtin -o z a.c -lm && ./z
cc -O3 -fno-builtin -o z a.c -lm && ./z

   ur ui: 21.588761 7.006300
     x y: 9.5623927 1.4993777
  csinhf: 5.07348328e+02 7.09178613e+03
dp_csinh: 5.07348986e+02 7.09178955e+03
   sinhf: 7.10991113e+03
    cosf: 7.13578984e-02

Max ULP: 23.061242
Count: 39          (# of ULP that exceeds 21)

Things are much worse than this toy program shows.
My test program used in development of libm is giving

Restrict x < 10

./testf -u -X 10
Max ULP Re: 136628.340239
Max ULP Im: 1891176.003955

Restrict c < 50
./testf -u -X 10
Max ULP Re: 3615923.332529
Max ULP Im: 13677733.591783

/*
 * Compute 1 million valus of csinhf() and then compute the ULP for
 * for the real and imaginary parts.
 */
#include <complex.h>
#include <float.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>

/* Return 0 <= x < 1. */
double
ranged(void)
{

   union {
      double x;
      struct {
         uint32_t lo;
         uint32_t hi;
      } u;
   } v;
   v.u.hi = (uint32_t)random();
   v.u.hi = ((v.u.hi << 11) >> 11) | 0x3ff00000;
   v.u.lo = (uint32_t)random();
   return (v.x - 1);
}

float
rangef(void)
{

        float s;
        s = (float)ranged();
        return (s);
}

/* Double precision csinh() without using C's double complex.s */
void
dp_csinh(double x, double y, double *re, double *im)
{
   double c, s;
   sincos(y, &s, &c);
   *re = sinh(x) * c;
   *im = cosh(x) * s;
}

/* ULP estimate. */
double
ulpfd(float app, double acc)
{
        int n;
        double f;
        f = frexp(acc, &n);
        f = fabs(acc - app);
        f = ldexp(f, FLT_MANT_DIG - n);
        return (f);     
}

int
main(void)
{
   double re, im, u, ur, ui;
   float complex f;
   float x, y;
   int cnt, i;

   srandom(19632019);

   ur = ui = 0;

   for (cnt = 0, i = 0; i < 10000000; i++) {
      x = rangef() + 9;
      y = rangef() + 0.5;
      f = csinhf(CMPLXF(x,y));
      dp_csinh((double)x, (double)y, &re, &im);
      ur = ulpfd(crealf(f), re);
      if (ur > u) u = ur;
      ui = ulpfd(cimagf(f), im);
      if (ui > u) u = ui;
      if (ur > 21 || ui > 21) {
         printf("   ur ui: %f %f\n", ur, ui);
         printf("     x y: %.7f %.7f\n", x, y);
         printf("  csinhf: %.8e %.8e\n", crealf(f), cimagf(f));
         printf("dp_csinh: %.8le %.8le\n", re, im);
         printf("   sinhf: %.8e\n", sinhf(x));
         printf("    cosf: %.8e\n\n", cosf(y));
         cnt++;
      }
   }
   printf("Max ULP: %f\n", u);
   printf("Count: %d\n", cnt);
   return (0);
}


-- 
Steve