Tobias Nygren <tnn%NetBSD.org@localhost> writes:
> On Thu, 11 Aug 2022 14:18:20 -0400
> Greg Troxel <gdt%lexort.com@localhost> wrote:
>
>> Has anyone else been having this problem?
>
> It's a known defect in NetBSD's libm.
I didn't know, but:
https://gnats.netbsd.org/52270
Thanks for the comments - a bunch of points I didn't realize became
clearer.
>> Any reason I shouldn't just define it and call log1p, not worrying about
>> those extra bits, to get around this (in my local build)?
>
> If you want to lose less precision I believe it is better
> to do "logl(1.0l + x)" than "(long double)log1p(x)".
Interesting idea, but I guess there's the missing bits vs the log1p
gain.
> But with homeassistant not being a scientific software I
> guess it don't matter much. Are you sure the issue is with
> numpy and not scipy? There is a patch for this issue in
> pkgsrc scipy.
Indeed I am 99.9% sure it isn't being called; I just need to have it
imported so that PyTurboJPEG can do some sort of jpeg/mpeg procesing. I
didn't actually miss the routine, so much as have the shlib use fail.
I'm pretty sure there is no scipy involved. The symbol log1pl appears
as U (via nm -g) in core/_multiarray_umath.so. I have the same numpy
version in virtualenv as is in pkgsrc and I can't explain why the pkgsrc
one succeeds at "import" time. But the virtualenv one imports ok, too,
and really there are multiple _foo.so there and it may be that they get
loaded later depending on what you do.
Fixing pkgsrc doesn't solve my issue; patches really need to be filed
upstream so that use outside of pkgsrc is ok, which is what happens in
pip virtualenv. But I can see that upstream might not want this, saying
they don't want patches
Why do we think it's useful/ok to declare the symbol in math.h if it's
not provided? That means you can't even write configure tests and avoid
using it. However it's probably better to fix it right.
>> Anyone up for doing it right?
>
> Does FreeBSD have an implementation in their libm?
That's where I was going to look...
> POSIX permits to just use the naive logl(1+x) implementation but
> the existence of log1p implies there exists a better precision
> numerical algorithm for x near 0.
Indeed, and I suppose we could have our compiler implement long double
as double and be legit, but at this point that's an ABI change.
I'm not suggesting committing, but this patch in my tree causes
homeassistant's use of from-pip numpy to work. This is a a terrible
hack and not cleaned up; on reviewing I se that the log2l and expm1l
cases are wrong! It's in b_exp.c as a random C file because the C file
containing log1p and log1pf is not compiled because the i387 .S version
is used instead.
But, we might consider this, cleaned up, or maybe log1l(1.0l+x), if we
can't steal something better from FreeBSD.
Index: lib/libm/src/b_exp.c
===================================================================
RCS file: /cvsroot/src/lib/libm/src/b_exp.c,v
retrieving revision 1.1
diff -u -p -r1.1 b_exp.c
--- lib/libm/src/b_exp.c 5 May 2012 17:54:14 -0000 1.1
+++ lib/libm/src/b_exp.c 12 Aug 2022 10:27:15 -0000
@@ -135,3 +135,41 @@ __exp__D(double x, double c)
/* exp(INF) is INF, exp(+big#) overflows to INF */
return( finite(x) ? scalb(1.0,5000) : x);
}
+
+
+/*
+ * Fake up long double version to let numpy build, and blow off it
+ * being more accurate. Put it here because s_log1p.c is really in
+ * asm.
+ * gdt, 20220811
+ */
+long double
+log1pl(long double x)
+{
+ long double y;
+
+ y = log1p(x);
+
+ return y;
+}
+
+long double
+log2l(long double x)
+{
+ long double y;
+
+ y = log2l(x);
+
+ return y;
+}
+
+long double
+expm1l(long double x)
+{
+ long double y;
+
+ y = expm1l(x);
+
+ return y;
+}
+
Attachment:
signature.asc
Description: PGP signature