<div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote">On Sat, Apr 26, 2014 at 6:37 PM, Matthew Brett <span dir="ltr"><<a href="mailto:matthew.brett@gmail.com" target="_blank">matthew.brett@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">Hi,<br>
<div><div><br>
On Wed, Apr 23, 2014 at 11:59 AM, Matthew Brett <<a href="mailto:matthew.brett@gmail.com" target="_blank">matthew.brett@gmail.com</a>> wrote:<br>
> Hi,<br>
><br>
> On Wed, Apr 23, 2014 at 1:43 AM, Nathaniel Smith <<a href="mailto:njs@pobox.com" target="_blank">njs@pobox.com</a>> wrote:<br>
>> On Wed, Apr 23, 2014 at 6:22 AM, Matthew Brett <<a href="mailto:matthew.brett@gmail.com" target="_blank">matthew.brett@gmail.com</a>> wrote:<br>
>>> Hi,<br>
>>><br>
>>> I'm exploring Mingw-w64 for numpy building, and I've found it gives a<br>
>>> slightly different answer for 'exp' than - say - gcc on OSX.<br>
>>><br>
>>> The difference is of the order of the eps value for the output number<br>
>>> (2 * eps for a result of ~2.0).<br>
>>><br>
>>> Is accuracy somewhere specified for C functions like exp? Or is<br>
>>> accuracy left as an implementation detail for the C library author?<br>
>><br>
>> C99 says (sec 5.2.4.2.2) that "The accuracy of the floating point<br>
>> operations ... and of the library functions in <math.h> and<br>
>> <complex.h> that return floating point results is implemenetation<br>
>> defined. The implementation may state that the accuracy is unknown."<br>
>> (This last sentence is basically saying that with regard to some<br>
>> higher up clauses that required all conforming implementations to<br>
>> document this stuff, saying "eh, who knows" counts as documenting it.<br>
>> Hooray for standards!)<br>
>><br>
>> Presumably the accuracy in this case is a function of the C library<br>
>> anyway, not the compiler?<br>
><br>
> Mingw-w64 implementation is in assembly:<br>
><br>
> <a href="http://sourceforge.net/p/mingw-w64/code/HEAD/tree/trunk/mingw-w64-crt/math/exp.def.h" target="_blank">http://sourceforge.net/p/mingw-w64/code/HEAD/tree/trunk/mingw-w64-crt/math/exp.def.h</a><br>
><br>
>> Numpy has its own implementations for a<br>
>> bunch of the math functions, and it's been unclear in the past whether<br>
>> numpy or the libc implementations were better in any particular case.<br>
><br>
> I only investigated this particular value, in which case it looked as<br>
> though the OSX value was closer to the exact value (via sympy.mpmath)<br>
> - by ~1 unit-at-the-last-place. This was causing a divergence in the<br>
> powell optimization path and therefore a single scipy test failure. I<br>
> haven't investigated further - was wondering what investigation I<br>
> should do, more than running the numpy / scipy test suites.<br>
<br>
</div></div>Investigating further, with this script:<br>
<br>
<a href="https://gist.github.com/matthew-brett/11301221" target="_blank">https://gist.github.com/matthew-brett/11301221</a><br>
<br>
The following are tests of np.exp accuracy for input values between 0<br>
and 10, for numpy 1.8.1.<br>
<br>
If np.exp(x) performs perfectly, it will return the nearest floating<br>
point value to the exact value of exp(x). If it does, this scores a<br>
zero for error in the tables below. If 'proportion of zeros' is 1 -<br>
then np.exp performs perfectly for all tested values of exp (as is the<br>
case for linux here).<br>
<br>
OSX 10.9<br>
<br>
Proportion of zeros: 0.99789<br>
Sum of error: 2.15021267458e-09<br>
Sum of squared error: 2.47149370032e-14<br>
Max / min error: 5.96046447754e-08 -2.98023223877e-08<br>
Sum of squared relative error: 5.22456992025e-30<br>
Max / min relative error: 2.19700100681e-16 -2.2098803255e-16<br>
eps: 2.22044604925e-16<br>
Proportion of relative err >= eps: 0.0<br>
<br>
Debian Jessie / Sid<br>
<br>
Proportion of zeros: 1.0<br>
Sum of error: 0.0<br>
Sum of squared error: 0.0<br>
Max / min error: 0.0 0.0<br>
Sum of squared relative error: 0.0<br>
Max / min relative error: 0.0 0.0<br>
eps: 2.22044604925e-16<br>
Proportion of relative err >= eps: 0.0<br>
<br>
Mingw-w64 Windows 7<br>
<br>
Proportion of zeros: 0.82089<br>
Sum of error: 8.08415331122e-07<br>
Sum of squared error: 2.90045099615e-12<br>
Max / min error: 5.96046447754e-08 -5.96046447754e-08<br>
Sum of squared relative error: 4.18466468175e-28<br>
Max / min relative error: 2.22041308226e-16 -2.22042100773e-16<br>
eps: 2.22044604925e-16<br>
Proportion of relative err >= eps: 0.0<br>
<br>
Take-home : exp implementation for mingw-w64 is exactly (floating<br>
point) correct 82% of the time, and one unit-at-the-last-place off for<br>
the rest [1]. OSX is off by 1 ULP only 0.2% of the time.<br></blockquote><div><br></div><div><br></div><div>Windows 64 with MKL</div><div><br></div><div><div>\WinPython-64bit-3.3.2.2\python-3.3.2.amd64>python "E:\Josef\eclipsegworkspace\statsmodels-git\local_scripts\local_scripts\try_exp_error.py"</div>
<div>Proportion of zeros: 0.99793</div><div>Sum of error: -2.10546855506e-07</div><div>Sum of squared error: 3.33304327526e-14</div><div>Max / min error: 5.96046447754e-08 -5.96046447754e-08</div><div>Sum of squared relative error: 4.98420694339e-30</div>
<div>Max / min relative error: 2.20881302691e-16 -2.18321571939e-16</div><div>eps: 2.22044604925e-16</div><div>Proportion of relative err >= eps: 0.0</div></div><div><br></div><div><br></div><div>Windows 32 bit python with official MingW binaries</div>
<div><br></div><div><div>Python 2.7.1 (r271:86832, Nov 27 2010, 18:30:46) [MSC v.1500 32 bit (Intel)] on win32</div><div><br></div><div>Proportion of zeros: 0.99464</div><div>Sum of error: -3.91621083118e-07</div><div>Sum of squared error: 9.2239247812e-14</div>
<div>Max / min error: 5.96046447754e-08 -5.96046447754e-08</div><div>Sum of squared relative error: 1.3334972729e-29</div><div>Max / min relative error: 2.21593462148e-16 -2.2098803255e-16</div><div>eps: 2.22044604925e-16</div>
<div>Proportion of relative err >= eps: 0.0</div></div><div><br></div><div> </div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<br>
Is mingw-w64 accurate enough? Do we have any policy on this?<br></blockquote><div><br></div><div><div>I wouldn't worry about a missing or an extra eps in our applications, but the competition is more accurate.</div>
<div><br></div><div>Josef</div></div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<br>
Cheers,<br>
<br>
Matthew<br>
<br>
[1] <a href="http://matthew-brett.github.io/pydagogue/floating_error.html" target="_blank">http://matthew-brett.github.io/pydagogue/floating_error.html</a><br>
<div><div>_______________________________________________<br>
NumPy-Discussion mailing list<br>
<a href="mailto:NumPy-Discussion@scipy.org" target="_blank">NumPy-Discussion@scipy.org</a><br>
<a href="http://mail.scipy.org/mailman/listinfo/numpy-discussion" target="_blank">http://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br>
</div></div></blockquote></div><br></div></div>