This is apropos numpy#6127. The fix is to make the functions floor_division
and remainder consistent, i.e.,

    b * floor_division(a, b) + remainder(a, b) == a

Previous to this fix remainder was computed a the C level using the '%'
operator, and the result was not always consistent with the floor
function. The current approach is to compute the remainder using

    b * (a/b - floor(a/b))

which is both consistent with the Python '%' operator and numerically
consistent with floor_division implemented using the floor function.

Closes numpy#6127.