[Numpy-discussion] Linking other libm-Implementation

Nils Becker nilsc.becker at gmail.com
Wed Feb 10 14:01:31 EST 2016


2016-02-09 18:02 GMT+01:00 Gregor Thalhammer <gregor.thalhammer at gmail.com>:
>> It is not suitable as a standard for numpy.
>
> Why should numpy not provide fast transcendental math functions? For
linear algebra it supports fast implementations, even non-free (MKL).
Wouldn’t it be nice if numpy outperforms C?

Floating point operations that make use of vector extensions of modern
processors may behave subtly different. This especially concerns floating
point exceptions, where sometimes silent infinities are generated instead
of raising a divide-by-zero exception (best description I could find on the
spot:
https://randomascii.wordpress.com/2012/04/21/exceptional-floating-point/,
also see the notes on C99-compliance of the new vector expressions in
glibc: https://sourceware.org/glibc/wiki/libmvec).
I think the default in numpy should strive to be mostly standard compliant.
But of course an option to activate vector math operations would be nice -
if that is necessary with packages like numexpr is another question.
One other point is the extended/long double type which is normally not
supported by those libraries (as vector extensions cannot handle them).

> Intel publishes accuracy/performance charts for VML/MKL:
>
https://software.intel.com/sites/products/documentation/doclib/mkl/vm/functions/_accuracyall.html
>
> For GNU libc it is more difficult to find similarly precise data, I only
could find:
>
http://www.gnu.org/software/libc/manual/html_node/Errors-in-Math-Functions.html

On Tue, Feb 9, 2016 at 7:06 AM, Daπid <davidmenhur at gmail.com> wrote:
> I did some digging, and I found this:
>
>
http://julia-programming-language.2336112.n4.nabble.com/Is-the-accuracy-of-Julia-s-elementary-functions-exp-sin-known-td32736.html

Thank you for looking that up! I did not knew about the stuff published by
Intel yet.

2016-02-09 20:13 GMT+01:00 Matthew Brett <matthew.brett at gmail.com>:
> So GNU libm has max error <= 0.5 ULP, openlibm has <= 1 ULP, and OSX
> is (almost always) somewhere in-between.
>
> So, is <= 1 ULP good enough?

Calculating transcendental functions correctly rounded is very, very hard
and to my knowledge there is no complete libm implementation that
guarantees the necessary accuracy for all possible inputs. One effort
was/is the Correctly Rounded LibM (crlibm [1]) which tried to prove the
accuracy of their algorithms. However, the performance impact to achieve
that last ulp in all rounding modes can be severe.
Assessing accuracy of a function implementation is hard. Testing all
possible inputs is not feasible (2^32/64 for single/double) and proving
accuracy bounds may be even harder.
Most of the time one samples accuracy with random numbers from a certain
range. This generates tables like the ones for GNU libm or Intel.
This is a kind of "faithful" accuracy as you believe that the accuracy you
tested on a sample extends to the whole argument range. The error in worst
case may be (much) bigger.

That being said, I believe the values given by GNU libm for example are
very trustworthy. libm is not always correctly rounded (which would be <=
0.5ulp in round-to-nearest), however, the error bounds given in the table
seem to cover all worst cases.
Common single-argument functions (sin, cos) are correctly rounded and even
complex two-argument functions (cpow) are at most 5ulp off. I do not think
that other implementations are more accurate.
So libm is definitely good enough, accuracy-wise.

In any case I would like to build a testing framework to compare some libms
and check accuracy/performance (at least Intel has a history of
underestimating their error bounds in transcendental functions [2]). crlibm
offers worst-case arguments for some functions which could be used to
complement randomized sampling. Maybe I have some time in the next weeks...

[1] http://lipforge.ens-lyon.fr/www/crlibm/
[2]
https://randomascii.wordpress.com/2014/10/09/intel-underestimates-error-bounds-by-1-3-quintillion/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20160210/e08292b7/attachment.html>


More information about the NumPy-Discussion mailing list