Date: Tue, 12 Mar 2019 19:45:41 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: freebsd-toolchain@freebsd.org, freebsd-current@freebsd.org Subject: Optimization bug with floating-point? Message-ID: <20190313024506.GA31746@troutmask.apl.washington.edu>
next in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20190313024506.GA31746>