[Numpy-discussion] Slightly off-topic - accuracy of C exp function?
Matthew Brett
matthew.brett at gmail.com
Sat Apr 26 18:37:25 EDT 2014
Hi,
On Wed, Apr 23, 2014 at 11:59 AM, Matthew Brett <matthew.brett at gmail.com> wrote:
> Hi,
>
> On Wed, Apr 23, 2014 at 1:43 AM, Nathaniel Smith <njs at pobox.com> wrote:
>> On Wed, Apr 23, 2014 at 6:22 AM, Matthew Brett <matthew.brett at gmail.com> wrote:
>>> Hi,
>>>
>>> I'm exploring Mingw-w64 for numpy building, and I've found it gives a
>>> slightly different answer for 'exp' than - say - gcc on OSX.
>>>
>>> The difference is of the order of the eps value for the output number
>>> (2 * eps for a result of ~2.0).
>>>
>>> Is accuracy somewhere specified for C functions like exp? Or is
>>> accuracy left as an implementation detail for the C library author?
>>
>> C99 says (sec 5.2.4.2.2) that "The accuracy of the floating point
>> operations ... and of the library functions in <math.h> and
>> <complex.h> that return floating point results is implemenetation
>> defined. The implementation may state that the accuracy is unknown."
>> (This last sentence is basically saying that with regard to some
>> higher up clauses that required all conforming implementations to
>> document this stuff, saying "eh, who knows" counts as documenting it.
>> Hooray for standards!)
>>
>> Presumably the accuracy in this case is a function of the C library
>> anyway, not the compiler?
>
> Mingw-w64 implementation is in assembly:
>
> http://sourceforge.net/p/mingw-w64/code/HEAD/tree/trunk/mingw-w64-crt/math/exp.def.h
>
>> Numpy has its own implementations for a
>> bunch of the math functions, and it's been unclear in the past whether
>> numpy or the libc implementations were better in any particular case.
>
> I only investigated this particular value, in which case it looked as
> though the OSX value was closer to the exact value (via sympy.mpmath)
> - by ~1 unit-at-the-last-place. This was causing a divergence in the
> powell optimization path and therefore a single scipy test failure. I
> haven't investigated further - was wondering what investigation I
> should do, more than running the numpy / scipy test suites.
Investigating further, with this script:
https://gist.github.com/matthew-brett/11301221
The following are tests of np.exp accuracy for input values between 0
and 10, for numpy 1.8.1.
If np.exp(x) performs perfectly, it will return the nearest floating
point value to the exact value of exp(x). If it does, this scores a
zero for error in the tables below. If 'proportion of zeros' is 1 -
then np.exp performs perfectly for all tested values of exp (as is the
case for linux here).
OSX 10.9
Proportion of zeros: 0.99789
Sum of error: 2.15021267458e-09
Sum of squared error: 2.47149370032e-14
Max / min error: 5.96046447754e-08 -2.98023223877e-08
Sum of squared relative error: 5.22456992025e-30
Max / min relative error: 2.19700100681e-16 -2.2098803255e-16
eps: 2.22044604925e-16
Proportion of relative err >= eps: 0.0
Debian Jessie / Sid
Proportion of zeros: 1.0
Sum of error: 0.0
Sum of squared error: 0.0
Max / min error: 0.0 0.0
Sum of squared relative error: 0.0
Max / min relative error: 0.0 0.0
eps: 2.22044604925e-16
Proportion of relative err >= eps: 0.0
Mingw-w64 Windows 7
Proportion of zeros: 0.82089
Sum of error: 8.08415331122e-07
Sum of squared error: 2.90045099615e-12
Max / min error: 5.96046447754e-08 -5.96046447754e-08
Sum of squared relative error: 4.18466468175e-28
Max / min relative error: 2.22041308226e-16 -2.22042100773e-16
eps: 2.22044604925e-16
Proportion of relative err >= eps: 0.0
Take-home : exp implementation for mingw-w64 is exactly (floating
point) correct 82% of the time, and one unit-at-the-last-place off for
the rest [1]. OSX is off by 1 ULP only 0.2% of the time.
Is mingw-w64 accurate enough? Do we have any policy on this?
Cheers,
Matthew
[1] http://matthew-brett.github.io/pydagogue/floating_error.html
More information about the NumPy-Discussion
mailing list