2016-02-09 18:02 GMT+01:00 Gregor Thalhammer <gregor.thalham...@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 <davidmen...@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.br...@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/
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org https://mail.scipy.org/mailman/listinfo/numpy-discussion