Re: [Mingw-w64-public] Wrong quotient results of `remquo()`?

[lhmouse]

Reliable results could usually be generated using GCC's constant folding.
For example, the following program:

#include 
int main(){
const double x = 10.001000, y = -1.299000;
int quo;
double rem = __builtin_remquo(x, y, );
printf("rem = %f, quo = %d\n", rem, quo);
}

after compiled with `gcc test.c -O3 -masm=intel -S`, produces the following 
assembly:

movsd   xmm0, QWORD PTR .LC0[rip]
lea rcx, .LC1[rip]
mov r8d, -8# <== folded constant goes here
movapd  xmm1, xmm0
movqrdx, xmm0
callprintf

And yes, the result -8 is correct. 0 isn't.

It doesn't matter whether we return -16 or -8 here, as both ISO C and POSIX
only require three bits at least. Here they are both conforming as long as
we document our `remquo()` as returning only three bits into `quo`.

--   
Best regards,
lh_mouse
2016-09-12

-
发件人："K. Frank" 
发送日期：2016-09-12 22:23
收件人：mingw64
抄送：
主题：Re: [Mingw-w64-public] Wrong quotient results of `remquo()`?

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  wrote:
> More likely a bug in mingw-w64:
>
> #include 
> #include 
> volatile double x = 10.001000, y = -1.299000;
> int main(){
> int quo;
> double rem = remquo(x, y, );
> 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 

Re: [Mingw-w64-public] Wrong quotient results of `remquo()`?

[K. Frank]

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  wrote:
> More likely a bug in mingw-w64:
>
> #include 
> #include 
> volatile double x = 10.001000, y = -1.299000;
> int main(){
> int quo;
> double rem = remquo(x, y, );
> 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