Hello Lefty!
I do think you have found a bug here, and it does appear to
be in the mingw-w64 code. Disclaimer: I don't understand
this completely.
Further comments in line, below.
On Tue, Sep 6, 2016 at 11:52 PM, lhmouse <lh_mo...@126.com> wrote:
> More likely a bug in mingw-w64:
>
> #include <stdio.h>
> #include <math.h>
> volatile double x = 10.001000, y = -1.299000;
> int main(){
> int quo;
> double rem = remquo(x, y, &quo);
> printf("rem = %f, quo = %d\n", rem, quo);
> }
>
> With mingw-w64 this program gives the following output:
>
> E:\Desktop>gcc test.c
>
> E:\Desktop>a
> rem = -0.391000, quo = 0
I get the same result as you in a c++ test program using:
g++ (x86_64-posix-seh-rev0, Built by MinGW-W64 project) 4.9.2
> However, according to ISO C11 draft:
>
>> 7.12.10.3 The remquo functions
>> 2 The remquo functions compute the same remainder as the remainder
>> functions. In
>> the object pointed to by quo they store a value whose sign is the sign of
>> x/y and whose
>> magnitude is congruent modulo 2n to the magnitude of the integral quotient
>> of x/y, where
>> n is an implementation-defined integer greater than or equal to 3.
>
> the value stored in `quo` must have the same sign with `x/y`.
>
> In the above example, since `x/y`, which is about -7.699, is negative,
> returning a non-negative value (zero in the above example) should be a bug.
I agree with lh_mouse's reading of the standard, and that
quo should be negative to match the sign of x / y.
Here is my imperfect analysis of what is going on.
First, I found a copy of remquo.S here:
https://sourceforge.net/p/mingw-w64/code/6570/tree//stable/v1.x/mingw-w64-crt/math/remquo.S
(but I don't understand it).
I also found a "softmath" copy of remquo.c here:
https://github.com/msys2/mingw-w64/blob/master/mingw-w64-crt/math/softmath/remquo.c
(I have no idea whether remquo.c is equivalent in detail to
remquo.S.)
#include "softmath_private.h"
double remquo(double x, double y, int *quo)
{
double r;
if (isnan(x) || isnan(y)) return NAN;
if (y == 0.0)
{
errno = EDOM;
__mingw_raise_matherr (_DOMAIN, "remquo", x, y, NAN);
return NAN;
}
r = remainder(x, y);
if (quo) *quo = (int)((x - r) / y) % 8;
return r;
}
First, the expression "(int)((x - r) / y)" is undefined
behavior when (x - r) / y is too large for an int. (This
can easily happen with floats and doubles.) (remquo.S
uses the intel floating-point instruction fprem1, and
therefore -- if written correctly -- should not have this
problem.)
But ignoring the possible integer overflow, the error here,
which is the result lh_mouse gets in his test, is that if
(int)((x - r) / y)
is a multiple of 8, then
(int)((x - r) / y) % 8
will evaluate to zero, losing the sign information.
In lh_mouse's test case
x / y = -7.699
which rounds-to-nearest to -8, which equals 0 mod 8.
How might one fix this (in c code)? My reading of the
standard says that quo doesn't have to be exactly the
last three bits -- or even the last n bits -- of
(int)((x - r) / y)
Rather, it only has to be congruent to this mod 8.
So (ignoring overflow in the integer conversion), one
could do something like this:
r = remainder(x, y);
if (quo) *quo = (int)((x - r) / y) % 8;
if (quo && *quo == 0 && x/y < 0.0) *quo = -16;
return r;
(Here, we are deeming the sign of 0 to be positive. I don't
know whether this would be language-lawyer consistent with
the standard, but for two's-complement integers, zero does
have the "positive" sign bit.)
The core problem is that the integer quotient could be a
large, negative multiple of 8 (and not even fit in an int).
In such a case, quo has to preserve the sign of x/y, so it
can't be 0. but there's no obvious favored (congruent mod 8)
approximation to x/y (other than 0, which we can't use), so
you might as well just use -16.
Should you also use +16 for quo when x/y is a positive
multiple of 8? My reading of the standard says that you
could, but that seems a little weird to me, so maybe 0
is more natural in the positive case.
A couple of cautions: I have no idea whether remquo.S
behaves the same way as remquo.c. However, remquo.c
does give the same incorrect value (0) for quo that
lh_mouse reports and that I see in my test. Also, I
assume that my test program ends up using remquo.S,
but I really haven't any idea how to check this.
In any event, remquo.c (at least the version I found
on line) is wrong on two counts -- it shows the bug
lh_mouse found, and it can have undefined behavior in
the integer conversion.
My guess is that lh_mouse's test and mine both end up
using remquo.S (Why would we end up using softmath?),
so presumably remquo.S has the quo == 0 bug (but, if I
had to guess, doesn't have the integer-overflow bug).
> Best regards,
> lh_mouse
> 2016-09-07
Best regards.
K. Frank
------------------------------------------------------------------------------
_______________________________________________
Mingw-w64-public mailing list
Mingw-w64-public@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public
