Possible miscompilation issue with floating point operations
[Please Cc me on replies, I'm not subscribed to the list]
Dear mips porters,
I attach a C source code that, when compiled with gcc-5, crashes with
"Illegal instruction" on mips, while it works correctly when compiled
with gcc-4.9 [tested on minkus.debian.org].
The problem looks like a miscompilation (the code only does some basic
floating-point operations), but it was pointed out by Julien Cristau
that it could also be a kernel issue (since floating point instructions
are emulated on Debian mips hardware).
What is your take on this? Should I file a bug against gcc?
Note that the attached program is a reduced testcase for #810357
(crashes of r-base with SIGILL on mips), that was solved by compiling
r-base with gcc-4.9. It may also be related, though this remains to be
verified, to the FTBFS of atlas on mips [0].
Thanks in advance for your feedback.
[0] https://buildd.debian.org/status/fetch.php?pkg=atlas&arch=mips&ver=3.10.2-9%2Bb1&stamp=1452968364
--
.''`. Sébastien Villemot
: :' : Debian Developer
`. `' http://sebastien.villemot.name
`- GPG Key: 4096R/381A7594
#include <math.h>
static void
machar(int *ibeta, int *it, int *irnd, int *ngrd, int *machep, int *negep,
int *iexp, int *minexp, int *maxexp, double *eps,
double *epsneg, double *xmin, double *xmax)
{
volatile double a, b, beta, betain, betah, one,
t, temp, tempa, temp1, two, y, z, zero;
int i, itemp, iz, j, k, mx, nxres;
one = 1;
two = one+one;
zero = one-one;
/* determine ibeta, beta ala malcolm. */
a = one;
do {
a = a + a;
temp = a + one;
temp1 = temp - a;
}
while(temp1 - one == zero);
b = one;
do {
b = b + b;
temp = a + b;
itemp = (int)(temp - a);
}
while (itemp == 0);
*ibeta = itemp;
beta = *ibeta;
/* determine it, irnd */
*it = 0;
b = one;
do {
*it = *it + 1;
b = b * beta;
temp = b + one;
temp1 = temp - b;
}
while(temp1 - one == zero);
*irnd = 0;
betah = beta / two;
temp = a + betah;
if (temp - a != zero)
*irnd = 1;
tempa = a + beta;
temp = tempa + betah;
if (*irnd == 0 && temp - tempa != zero)
*irnd = 2;
/* determine negep, epsneg */
*negep = *it + 3;
betain = one / beta;
a = one;
for(i=1 ; i<=*negep ; i++)
a = a * betain;
b = a;
for(;;) {
temp = one - a;
if (temp - one != zero)
break;
a = a * beta;
*negep = *negep - 1;
}
*negep = -*negep;
*epsneg = a;
if (*ibeta != 2 && *irnd != 0) {
a = (a * (one + a)) / two;
temp = one - a;
if (temp - one != zero)
*epsneg = a;
}
/* determine machep, eps */
*machep = -*it - 3;
a = b;
for(;;) {
temp = one + a;
if (temp - one != zero)
break;
a = a * beta;
*machep = *machep + 1;
}
*eps = a;
temp = tempa + beta * (one + *eps);
if (*ibeta != 2 && *irnd != 0) {
a = (a * (one + a)) / two;
temp = one + a;
if (temp - one != zero)
*eps = a;
}
/* determine ngrd */
*ngrd = 0;
temp = one + *eps;
if (*irnd == 0 && temp * one - one != zero)
*ngrd = 1;
/* determine iexp, minexp, xmin */
/* loop to determine largest i and k = 2**i such that */
/* (1/beta) ** (2**(i)) */
/* does not underflow. */
/* exit from loop is signaled by an underflow. */
i = 0;
k = 1;
z = betain;
t = one + *eps;
nxres = 0;
for(;;) {
y = z;
z = y * y;
/* check for underflow here */
a = z * one;
temp = z * t;
if (a+a == zero || fabs(z) >= y)
break;
temp1 = temp * betain;
if (temp1 * beta == z)
break;
i = i+1;
k = k+k;
}
if (*ibeta != 10) {
*iexp = i + 1;
mx = k + k;
}
else {
/* this segment is for decimal machines only */
*iexp = 2;
iz = *ibeta;
while (k >= iz) {
iz = iz * *ibeta;
iexp = iexp + 1;
}
mx = iz + iz - 1;
}
do {
/* loop to determine minexp, xmin */
/* exit from loop is signaled by an underflow */
*xmin = y;
y = y * betain;
/* check for underflow here */
a = y * one;
temp = y * t;
if (a+a == zero || fabs(y) >= *xmin)
goto L10;
k = k + 1;
temp1 = temp * betain;
}
while(temp1 * beta != y);
nxres = 3;
*xmin = y;
L10: *minexp = -k;
/* determine maxexp, xmax */
if (mx <= k + k - 3 && *ibeta != 10) {
mx = mx + mx;
*iexp = *iexp + 1;
}
*maxexp = mx + *minexp;
/* adjust irnd to reflect partial underflow */
*irnd = *irnd + nxres;
/* adjust for ieee-style machines */
if (*irnd == 2 || *irnd == 5)
*maxexp = *maxexp - 2;
/* adjust for non-ieee machines with partial underflow */
if (*irnd == 3 || *irnd == 4)
*maxexp = *maxexp - *it;
/* adjust for machines with implicit leading bit in binary */
/* significand, and machines with radix point at extreme */
/* right of significand. */
i = *maxexp + *minexp;
if (*ibeta == 2 && i == 0)
*maxexp = *maxexp - 1;
if (i > 20)
*maxexp = *maxexp - 1;
if (a != y)
*maxexp = *maxexp - 2;
*xmax = one - *epsneg;
if (*xmax * one != *xmax)
*xmax = one - beta * *epsneg;
*xmax = *xmax / (beta * beta * beta * *xmin);
i = *maxexp + *minexp + 3;
if (i>0)
for(j=1 ; j<=i ; j++) {
if (*ibeta == 2)
*xmax = *xmax + *xmax;
if (*ibeta != 2)
*xmax = *xmax * beta;
}
}
int
main(int argc, char **argv)
{
int ibeta, it, irnd, ngrd, machep, negep,
iexp, minexp, maxexp;
double eps, epsneg, xmin, xmax;
machar(&ibeta, &it, &irnd, &ngrd, &machep, &negep,
&iexp, &minexp, &maxexp, &eps,
&epsneg, &xmin, &xmax);
return 0;
}
Reply to: