Re: [Mingw-w64-public] Wrong quotient results of `remquo()`?
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, &quo); 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, &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 ca
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, &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 -