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 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. Unfortunately, I did not follow the (many different) bug reports on this issue originally and am confused why there was a need to change the implementation in the first place. numpy's "%" operator seems to work quite well on my system. Therefore, this mail may be rather unproductive as I am missing some information. Concerning your original question: Many elementary functions loose their mathematical properties when they are calculated correctly-rounded in floating point numbers [2]. We do not fix this for other functions, I would not fix it here. Cheers Nils [1] https://github.com/JuliaLang/openlibm/blob/master/src/e_fmod.c [2] np.exp(np.nextafter(1.0, 0.0)) < np.e -> False (Monotonicity lost in exp(x)).
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org https://mail.scipy.org/mailman/listinfo/numpy-discussion
