Weird result of floor division in numpy

Question

I stumbled upon a floor division result of a np.float32 or np.float64 , which I do not understand

I'm using numpy 1.15.4 in Python 3.6.7

>>> import numpy as np
>>> np.float32(0)//1
-0.0

I know, that i can workaround by abs() and such, but why am I getting a "-0.0" instead of "0.0" in the first place??

Answer 1

I suspect that numpy uses the divmod function to compute the floor division and the line causing this is here:

/* if div is zero ensure correct sign */
floordiv = (a / b > 0) ? 0.0@c@ : -0.0@c@;

Example:

>>> a = np.zeros(1)
>>> b = 1
>>> np.where(a/b > 0, 0.0, -0.0)
array([-0.])

Python's divmod function seems to handle this correctly, so they must be using different algorithms:

>>> divmod(0.0,1)
(0.0, 0.0)
>>> divmod(-0.0,1)
(-0.0, 0.0)

I looked into this a bit more and here is how python's divmod works for floating point numbers when div is zero (link):

/* div is zero - get the same sign as the true quotient */
floordiv = copysign(0.0, vx / wx); /* zero w/ sign of vx/wx */

and copysign() is defined as:

double
copysign(double x, double y)
{
    /* use atan2 to distinguish -0. from 0. */
    if (y > 0. || (y == 0. && atan2(y, -1.) > 0.)) {
        return fabs(x);
    } else {
        return -fabs(x);
    }
}

so the reason python is able to get this right and numpy isn't is that python uses atan2() to distinguish between -0.0 and +0.0.

Update: This issue will be fixed in the 1.17.0 release of numpy. You can see the release notes here.