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

2016-09-12 Thread 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" <kfrank2...@gmail.com>
发送日期: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 <lh_mo...@126.com> 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
t

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

2016-09-12 Thread 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


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

2016-09-06 Thread lhmouse
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

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.

--   
Best regards,
lh_mouse
2016-09-07




--
___
Mingw-w64-public mailing list
Mingw-w64-public@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public


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

2016-09-05 Thread lhmouse
The disagreement of glibc and mingw-w64 (in my opinion) is definitely glibc's 
bug:

lh_mouse@lhmouse-dev:~$ cat test3.c 
#include 
#include 

int main(){
double x = 10.001000;
double y = 0.701000;
int quo;
double rem = remquo(x, y, );
printf("%f %f %d %f\n", x, y, quo, rem);
}
lh_mouse@lhmouse-dev:~$ gcc test3.c -lm -O0 && ./a.out # use glibc
10.001000 0.701000 8 4.393000
lh_mouse@lhmouse-dev:~$ gcc test3.c -lm -O2 && ./a.out  # performs constant 
folding
10.001000 0.701000 14 0.187000
lh_mouse@lhmouse-dev:~$ 

The remainder of `remquo` from mingw-w64 seems all right.
However the value (or rather, the 3 least significant bits) returned in the 
third parameter
still seems problematic.

--   
Best regards,
lh_mouse
2016-09-06

-
发件人:"lhmouse"<lh_mo...@126.com>
发送日期:2016-09-05 23:08
收件人:mingw-w64-public,lhmouse
抄送:
主题:Re:  [Mingw-w64-public] Wrong quotient results of `remquo()`?

Found an example on cppreference:
http://en.cppreference.com/w/cpp/numeric/math/remquo

The example shows that, since `cos()` is periodic,
adding 1 * PI to its parameter doesn't change the result.

But, we can also say that, subtracting 1 * PI
from its parameter should not change the result either. However,
with mingw-w64 and MSVCRT, it DOES change the result,
as shown on the last line:

E:\Desktop>g++ test.cpp -std=c++14

E:\Desktop>a.exe
cos(pi * -0.25) = 0.707107
cos(pi * -1.25) = -0.707107
cos(pi * -1.25) = 0.707123
cos(pi * -10001.25) = -0.707117
cos(pi * -1.25) = 0.707107
cos(pi * -10001.25) = 0.707107

This could be a potential bug.

--   
Best regards,
lh_mouse
2016-09-05

-
发件人:"lhmouse"<lh_mo...@126.com>
发送日期:2016-09-05 22:27
收件人:mingw-w64-public
抄送:
主题:[Mingw-w64-public] Wrong quotient results of `remquo()`?

Hello guys,
I am testing my `remquo()` implementation when I find that `remquo`
on Linux (using glibc) and on Windows (using mingw-w64) generate
different results. I don't think this is the correct behavior. Any ideas?

The testcases in file `remquo.txt` the attached zip file was generated
on my VPS running Debian. MinGW-w64 is failing some of them:

E:\Desktop\remquo_test>gcc test.c -std=c99 && a.exe > nul
passed: 37864
failed: 2537

--
Best regards,
lh_mouse
2016-09-05
--

___
Mingw-w64-public mailing list
Mingw-w64-public@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public





--
___
Mingw-w64-public mailing list
Mingw-w64-public@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public


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

2016-09-05 Thread lhmouse
Found an example on cppreference:
http://en.cppreference.com/w/cpp/numeric/math/remquo

The example shows that, since `cos()` is periodic,
adding 1 * PI to its parameter doesn't change the result.

But, we can also say that, subtracting 1 * PI
from its parameter should not change the result either. However,
with mingw-w64 and MSVCRT, it DOES change the result,
as shown on the last line:

E:\Desktop>g++ test.cpp -std=c++14

E:\Desktop>a.exe
cos(pi * -0.25) = 0.707107
cos(pi * -1.25) = -0.707107
cos(pi * -1.25) = 0.707123
cos(pi * -10001.25) = -0.707117
cos(pi * -1.25) = 0.707107
cos(pi * -10001.25) = 0.707107

This could be a potential bug.

--   
Best regards,
lh_mouse
2016-09-05

-
发件人:"lhmouse"<lh_mo...@126.com>
发送日期:2016-09-05 22:27
收件人:mingw-w64-public
抄送:
主题:[Mingw-w64-public] Wrong quotient results of `remquo()`?

Hello guys,
I am testing my `remquo()` implementation when I find that `remquo`
on Linux (using glibc) and on Windows (using mingw-w64) generate
different results. I don't think this is the correct behavior. Any ideas?

The testcases in file `remquo.txt` the attached zip file was generated
on my VPS running Debian. MinGW-w64 is failing some of them:

E:\Desktop\remquo_test>gcc test.c -std=c99 && a.exe > nul
passed: 37864
failed: 2537

--
Best regards,
lh_mouse
2016-09-05
--

___
Mingw-w64-public mailing list
Mingw-w64-public@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public




--
___
Mingw-w64-public mailing list
Mingw-w64-public@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public


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

2016-09-05 Thread lhmouse
Hello guys,
I am testing my `remquo()` implementation when I find that `remquo`
on Linux (using glibc) and on Windows (using mingw-w64) generate
different results. I don't think this is the correct behavior. Any ideas?

The testcases in file `remquo.txt` the attached zip file was generated
on my VPS running Debian. MinGW-w64 is failing some of them:

E:\Desktop\remquo_test>gcc test.c -std=c99 && a.exe > nul
passed: 37864
failed: 2537

--
Best regards,
lh_mouse
2016-09-05
--
___
Mingw-w64-public mailing list
Mingw-w64-public@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public