> Date: Sat, 6 Sep 2014 14:15:32 -0400
> From: Daniel Dickman <didick...@gmail.com>
>
> according to the numpy developers asinhl on sparc64 might be buggy. I
> haven't worked out a test case yet but just reporting in case anyone else
> wants to take a look as well.
>
> bug report:
> https://github.com/numpy/numpy/issues/5026#issuecomment-54711361
Our asinhl(3) implementation is actually pretty much identical to the
one in glibc. The problem is that our sqrtl(3) implementation is
buggy, and this function is used to calculate asinhl(3) in certain
regimes. Fix forthcoming and with that fix, the testcase in the
bugreport no longer fails.
The fix below is from FreeBSD. I can't say I completely understand
it, but it makes some sense at least. It's a bit different from the
FreeBSD commit as they put in some unrelated "optimizations" and
departed from the existing coding style a bit.
ok?
Index: libc/arch/sparc64/fpu/fpu_sqrt.c
===================================================================
RCS file: /cvs/src/lib/libc/arch/sparc64/fpu/fpu_sqrt.c,v
retrieving revision 1.3
diff -u -p -r1.3 fpu_sqrt.c
--- libc/arch/sparc64/fpu/fpu_sqrt.c 26 Nov 2013 20:33:07 -0000 1.3
+++ libc/arch/sparc64/fpu/fpu_sqrt.c 10 Sep 2014 19:47:38 -0000
@@ -372,12 +372,12 @@ __fpu_sqrt(fe)
FPU_SUBCS(d2, x2, t2);
FPU_SUBCS(d1, x1, t1);
FPU_SUBC(d0, x0, t0);
- ODD_DOUBLE;
if ((int)d0 >= 0) {
- x0 = d0, x1 = d1, x2 = d2;
+ x0 = d0, x1 = d1, x2 = d2, x3 = d3;
q |= bit;
y2 |= 1;
}
+ ODD_DOUBLE;
while ((bit >>= 1) != 0) {
EVEN_DOUBLE;
t3 = y3 | bit;
@@ -386,7 +386,7 @@ __fpu_sqrt(fe)
FPU_SUBCS(d1, x1, t1);
FPU_SUBC(d0, x0, t0);
if ((int)d0 >= 0) {
- x0 = d0, x1 = d1, x2 = d2;
+ x0 = d0, x1 = d1, x2 = d2, x3 = d3;
q |= bit;
y3 |= bit << 1;
}
