[Numpy-discussion] Modulus (remainder) function corner cases

[Charles R Harris]

On Sun, Feb 14, 2016 at 12:35 PM, Nils Becker <nilsc.becker at gmail.com>
wrote:

> 2016-02-13 17:42 GMT+01:00 Charles R Harris <charlesr.harris at gmail.com>:
>
>> The Fortran modulo function, which is the same basic function as in my
>>> branch, does not specify any bounds on the result for floating numbers, but
>>> gives only the formula,  modulus(a, b) = a - b*floor(a/b), which has
>>> the advantage of being simple and well defined ;)
>>>
>>
>>
> In the light of the libm-discussion I spent some time looking at floating
> point functions and their accuracy. I would vote in favor of keeping an
> implementation that uses the fmod-function of the system library and bends
> it to adhere to the python convention (sign of divisor). There is probably
> a reason why the fmod-implementation is not as simple as "a - b*floor(a/b)"
> [1].
>
> One obvious problem with the simple expression arises when a/b = 0.0 in
> floating point. E.g.
>
> In [43]: np.__version__
> Out[43]: '1.10.4'
> In [44]: x = np.float64(1e-320)
> In [45]: y = np.float64(-1e10)
> In [46]: x % y # this uses libm's fmod on my system
> Out[46]: -10000000000.0 # == y, correctly rounded result in
> round-to-nearest
> In [47]: x - y*np.floor(x/y) # this here is the naive expression
> Out[47]: 9.9998886718268301e-321  # == x, wrong sign
>

But more accurate ;) Currently, this is actually clipped.

In [3]: remainder(x,y)
Out[3]: -0.0

In [4]: x - y*floor(x/y)
Out[4]: 9.9998886718268301e-321

In [5]: fmod(x,y)
Out[5]: 9.9998886718268301e-321



> There are other problems (a/b = inf in floating point). As I do not
> understand the implementation of fmod (for example in openlibm) in detail I
> cannot give a full list of corner cases.
>

?

<snip>

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20160214/7d7d5402/attachment.html>