> Date: Wed, 10 Sep 2014 21:52:30 +0200 (CEST)
> From: Mark Kettenis <[email protected]>
>
> > Date: Sat, 6 Sep 2014 14:15:32 -0400
> > From: Daniel Dickman <[email protected]>
> >
> > 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?
Now the fix for _Qp_sqrt() soft-float routine is enough to fix
asinhl(3), but only because gcc decides to directly invoke it inside
asinhl(3). But in other contexts it decides to actually emit a
sqrtl(3) library call, which invokes a totally different
implementation. And that other implementation, which lives in
libm/src/e_sqrtl.c, is buggy too.
So here is a diff that overrides this file for sparc64 with some code
that simply invokes _Qp_sqrtl(3).
ok?
Index: libm/Makefile
===================================================================
RCS file: /cvs/src/lib/libm/Makefile,v
retrieving revision 1.106
diff -u -p -r1.106 Makefile
--- libm/Makefile 18 Mar 2014 22:36:30 -0000 1.106
+++ libm/Makefile 10 Sep 2014 19:56:16 -0000
@@ -67,6 +67,7 @@ ARCH_SRCS = e_sqrt.c e_sqrtf.c s_fabsf.c
.PATH: ${.CURDIR}/arch/sparc
.elif (${MACHINE_ARCH} == "sparc64")
.PATH: ${.CURDIR}/arch/sparc64
+ARCH_SRCS = e_sqrtl.c
.elif (${MACHINE_ARCH} == "vax")
.PATH: ${.CURDIR}/arch/vax
NOIEEE_ARCH = n_argred.S n_infnan.S n_sqrt.S
Index: libm/arch/sparc64/e_sqrtl.c
===================================================================
RCS file: libm/arch/sparc64/e_sqrtl.c
diff -N libm/arch/sparc64/e_sqrtl.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ libm/arch/sparc64/e_sqrtl.c 10 Sep 2014 19:56:16 -0000
@@ -0,0 +1,29 @@
+/* $OpenBSD$ */
+/*
+ * Copyright (c) 2014 Mark Kettenis
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include <math.h>
+
+extern void _Qp_sqrt(long double *, long double *);
+
+long double
+sqrtl(long double x)
+{
+ long double y;
+
+ _Qp_sqrt(&y, &x);
+ return y;
+}
