On Mon, Oct 09, 2017 at 10:39:47PM +0200, Juan Francisco Cantero Hurtado wrote:
> Marc Feeley (Gambit Scheme) has been helping me with a bug on Gambit on
> Loongson. Apparently the bug is on our side.
>
> I've created this little test based on his code:
>
> #include <stdio.h>
> #include <math.h>
>
> int main()
> {
> double x = 0.5;
> double y = 1074.0;
> printf("x=%.20g y=%.20g pow(x,y)=%.20g\n",x,y,pow(x,y));
> return 0;
> }
>
> On OpenBSD/amd64, Linux/amd64, Linux/aarch64 and Linux/mips64:
> x=0.5 y=1074 pow(x,y)=4.9406564584124654418e-324
>
> On OpenBSD/loongson:
> x=0.5 y=1074 pow(x,y)=-4.9406564584124654418e-324
>
> Is the negative result expected?
No. On mips64, ldexp(3) seems to use a garbage value when figuring out
the sign of a denormal value. This applies to infinity values as well.
The patch below should fix the bug. The t3 register contains the binary
representation of the floating-point `x' parameter. I can also make
a regression test for the patch.
OK?
Index: arch/mips64/gen/ldexp.S
===================================================================
RCS file: src/lib/libc/arch/mips64/gen/ldexp.S,v
retrieving revision 1.7
diff -u -p -r1.7 ldexp.S
--- arch/mips64/gen/ldexp.S 27 Oct 2015 05:54:49 -0000 1.7
+++ arch/mips64/gen/ldexp.S 10 Oct 2017 11:52:56 -0000
@@ -143,20 +143,20 @@ LEAF(ldexp, 0)
xori t2, 1
1:
dmtc1 t2, $f0 # save denormalized result (LSW)
- bge v1, zero, 1f # should result be negative?
+ bge t3, zero, 1f # should result be negative?
neg.d $f0, $f0 # negate result
1:
j ra
7:
dmtc1 zero, $f0 # result is zero
- beq t0, zero, 1f # is result positive?
+ bge t3, zero, 1f # is result positive?
neg.d $f0, $f0 # negate result
1:
j ra
8:
dli t1, 0x7ff0000000000000 # result is infinity (MSW)
dmtc1 t1, $f0
- bge v1, zero, 1f # should result be negative infinity?
+ bge t3, zero, 1f # should result be negative infinity?
neg.d $f0, $f0 # result is negative infinity
1:
add.d $f0, $f0 # cause overflow faults if enabled
