pow() returns a negative result on loongson

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

pow() returns a negative result on loongson

Juan Francisco Cantero Hurtado
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?


--
Juan Francisco Cantero Hurtado http://juanfra.info

Reply | Threaded
Open this post in threaded view
|

Re: pow() returns a negative result on loongson

Paul Irofti-4
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, even though what you are computing is ridiculous :)

My first thought would on Loongson would be softfloat bug.
But I am not 100% sure it (still) uses softfloat.

Reply | Threaded
Open this post in threaded view
|

Re: pow() returns a negative result on loongson

Miod Vallat
> My first thought would on Loongson would be softfloat bug.
> But I am not 100% sure it (still) uses softfloat.

Loongson has never used softfloat. But all mips ports embed the
softfloat code in order to perform tricky computations the FPU won't
perform (mostly edge cases involving denormals and infinities).

So in this case, either the computation gets entirely done in hardware,
but the hardware yields a wrong results, and there won't likely be a fix
for this, or it asks the kernel for assistance, and there might be a
sign bug in there.

Reply | Threaded
Open this post in threaded view
|

Re: pow() returns a negative result on loongson

Visa Hankala-2
In reply to this post by Juan Francisco Cantero Hurtado
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

Reply | Threaded
Open this post in threaded view
|

Re: pow() returns a negative result on loongson

Juan Francisco Cantero Hurtado
On Tue, Oct 10, 2017 at 12:08:06PM +0000, Visa Hankala wrote:

> 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?

Works for me. Thanks for the quick fix.

>
> 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
>

--
Juan Francisco Cantero Hurtado http://juanfra.info