A comparison of the various approaches:

>>> a, b = 3, 2.5
>>> np.remainder(a, b)
0.49999999999999989
>>> b * (np.true_divide(a, b) - np.floor_divide(a, b))
0.49999999999999989
>>> a - (b * np.floor_divide(a, b))
0.5

I don't think 1 ulp counts as in issue in floating point. The new approach works better in general than the old.

I was a bit surprised here because the computation can be done exactly. But, well, agreed that it's not a major issue.
(for the record, CPython uses the platform fmod)

I think the problem is somewhat philosophical. You know that you meant 5 and 2.5 exactly, but the computer only knows that 5 and 2.5 are approximations of the true numbers, which may not be exactly represented. OTOH, I won't argue that further improvements can't be made. The original problem was to keep the results of floor division and remainder compatible, and when fmod (%) was used for the latter it was possible to have large inconsistencies. The current implementation guarantees (I think) that the remainder is always less in absolute value than the divisor and has the correct sign.

the computer only knows that 5 and 2.5 are approximations of the true numbers, which may not be exactly represented.

As @pitrou has mentioned, those particular numbers are exactly represented; a rational number with a small enough numerator and with a denominator that is a small enough power of 2 is represented exactly in floating point.

@argriffing You missed the point. All floating point numbers exactly represent themselves.

For fun, also

In [27]: 2.5 + 0.49999999999999989
Out[27]: 3.0

In [28]: 2.5 + 0.5
Out[28]: 3.0

I agree that an insignificant error in np.remainder is a small price for having fixed the coordination issues among functions like divmod/div/truediv/floordiv/remainder/etc.

While a small error is definitely preferable over the wrong answers given before, I can see it being disconcerting that

In [11]: np.remainder(6., 5.)
Out[11]: 0.99999999999999978

I'm slightly puzzled by the inconsistency with % and np.fmod:

In [12]: 6. % 5.
Out[12]: 1.0

In [15]: np.fmod(6., 5.)
Out[15]: 1.0

But in the end, it may just be a question of "pick your poison":

In [56]: a = 6.

In [57]: b = 5.

In [58]: b * (a / b - (a // b)) - a % b
Out[58]: -2.220446049250313e-16

In [59]: a = np.array(6.)

In [60]: b = np.array(5.)

In [61]: b * (a / b - (a // b)) - a % b
Out[61]: 0.0

As I was getting curious about what python itself does, I went ahead and looked for it -- not that easy to trace, but here it is: http://sources.debian.net/src/python2.7/2.7.11-3/Objects/floatobject.c/?hl=2083#L749

In the end, I think the premise is slightly different: python takes fmod as "god-given" and adjusts the floor division if necessary, while @charris takes floor division as "god-given" and uses it to calculate the remainder.

But it may be possible to remove the rounding problems, at the cost of an extra multiplication. Right now, we have (https://github.com/numpy/numpy/blob/master/numpy/core/src/umath/loops.c.src#L1710):

*((@type@ *)op1) = in2*(div - npy_floor@c@(div));

I think if we were to change this to

*((@type@ *)op1) = in2*div - in2*npy_floor@c@(div));

it would remove the "regression". At least

In [107]: a = np.array(6.)

In [108]: b = np.array(5.)

In [109]: b * (a / b - a // b)
Out[109]: 0.99999999999999978

In [110]: b * (a / b) - b * (a // b)
Out[110]: 1.0

(I checked on the command line, but not with actual change to C code, that this did not immediately bring back the original problem.)

Interesting, that may be better. What I liked about the original is that the abs value of the result is guaranteed to be less than abs(b) (I think). That may also be true of your new form (not sure). I think the original gives the smallest floating solution of (rem + b*floor(a/b)) = a, 0 <= sign(b)*rem < abs(b).

Hmm,

the IEEE standard also specifies that square root, remainder, and conversion between integer and floating-point be correctly rounded.

@mhvk I think that works.

I'm slightly puzzled by the inconsistency with % and np.fmod

Turns out that they have different definitions in terms of how they are computed.

Well, turns out my computation doesn't work for negative numbers close to zero. The problem is that ieee floating point precision is symmetric around zero, but floor is not. Basically, the python % is not numerically good, whereas the more common computer language versions based on (symmetric) truncation rather than floor are. OTOH, for many things the floor version is more desirable. The problem is less severe at larger distance from zero as the asymmetry is less pronounced. Hmm...

@charris: if we're having doubts at the 11th hour, should we pull it back to master and wait for 1.12?

Don't forget that in 1.10 one gets answers that are wrong, not just slightly off!

@njsmith It's an interesting problem ;) The tiny negative number bit is inherent, so present in the python version also. For exactly represented numbers using few bits, the ideal computation is a - b*floor(a/b), as floor should be correct and then all the rest of the arithmetic also. The problem I had with that is that sometimes the result would not have the right sign or be too big due to roundoff even though the errors were tiny. So all of these are fixable.