[Numpy-discussion] Test error with ATLAS, Windows 64 bit

[Matthew Brett]

On Mon, Apr 14, 2014 at 5:01 PM, Matthew Brett <matthew.brett at gmail.com> wrote:
> Hi,
>
> On Mon, Apr 14, 2014 at 3:55 PM, Charles R Harris
> <charlesr.harris at gmail.com> wrote:
>>
>>
>>
>> On Mon, Apr 14, 2014 at 4:40 PM, Matthew Brett <matthew.brett at gmail.com>
>> wrote:
>>>
>>> Hi,
>>>
>>> On Mon, Apr 14, 2014 at 2:58 PM, Charles R Harris
>>> <charlesr.harris at gmail.com> wrote:
>>> >
>>> >
>>> >
>>> > On Mon, Apr 14, 2014 at 3:38 PM, Matthew Brett <matthew.brett at gmail.com>
>>> > wrote:
>>> >>
>>> >> Hi,
>>> >>
>>> >> On Mon, Apr 14, 2014 at 2:02 PM, Julian Taylor
>>> >> <jtaylor.debian at googlemail.com> wrote:
>>> >> > The official numpy mingw binaries do not have all these math issues.
>>> >> > Only the VC builds do.
>>> >> > As mingw is fine the functions must be somewhere in the windows API
>>> >> > but
>>> >> > no-one has contributed a fix for the VC builds to numpy yet.
>>> >>
>>> >> I'm building with mingw-w64.
>>> >>
>>> >> It looks like this works as expected from this test:
>>> >>
>>> >> #include <math.h>
>>> >> #include <stdio.h>
>>> >>
>>> >> int main() {
>>> >>     double z;
>>> >>     z = expm1(-0.0);
>>> >>     printf("Result %f\n", z);
>>> >> }
>>> >>
>>> >> (prints -0).
>>> >>
>>> >> as does this (modified from core/src/npymath/npy_math.c.src):
>>> >>
>>> >> #include <stdio.h>
>>> >>
>>> >> double npy_expm1(double x)
>>> >> {
>>> >>     if (isinf(x) && x > 0) {
>>> >>         return x;
>>> >>     }
>>> >>     else {
>>> >>         const double u = exp(x);
>>> >>
>>> >>         if (u == 1.0) {
>>> >>             return x;
>>> >>         } else if (u - 1.0 == -1.0) {
>>> >>             return -1;
>>> >>         } else {
>>> >>             return (u - 1.0) * x/log(u);
>>> >>         }
>>> >>     }
>>> >> }
>>> >>
>>> >> int main() {
>>> >>     double z;
>>> >>     z = npy_expm1(-0.0);
>>> >>     printf("Result %f\n", z);
>>> >> }
>>> >>
>>> >> Sorry for my ignorance, but where does the `HAVE_EXPM1` symbol come
>>> >> from?
>>> >>
>>> >
>>> > Remember all configuration output at the beginning of the build? One of
>>> > those tries to link 'expm1{f, , l}' and sets the value of
>>> > 'HAVE_EXPM1{F,,L}'
>>> > accordingly. You can find the result in the `config.h` file in the build
>>> > directory.
>>>
>>> Ah - thanks - that got me going - somewhere.
>>>
>>> I noticed that HAVE_EXPM1 is not defined in any of OSX, my Windows
>>> build, Debian config.h - is that expected?
>>
>>
>> It's defined for me,  gcc version 4.8.2 on fedora 20. I'm surprised that
>> Debian doesn't show it, what compiler version does Debian provide?
>
> Yes, sorry - I had an old version of numpy on Debian.  It is defined for 1.8.1
>
> It's not defined on OSX because it is defined in pyconfig.h (via Python.h).
>
> I was wrong, it is defined in config.h for Windows, I'm not sure how I
> missed that.  Now I'm trying to understand why it's giving different
> results from the tiny test function above,

It looks as though mingw-w64 is at fault, and I was confused (still
am) because of the different behavior with double and a constant:

#include <stdio.h>
#include <math.h>

int main() {
    double z, i = -0.0;
    printf("With double %f=%f, with constant=%f\n",
   i, expm1(i), expm1(-0.));
}

gives:

With double -0.000000=0.000000, with constant=-0.000000

That was ugly to track down.

What is the right way to work round this (using the numpy version
instead of the system version I suppose)?

Cheers,

Matthew