On Sun, Feb 14, 2016 at 12:35 PM, Nils Becker <nilsc.bec...@gmail.com> wrote:
> 2016-02-13 17:42 GMT+01:00 Charles R Harris <charlesr.har...@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
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org https://mail.scipy.org/mailman/listinfo/numpy-discussion
