Modulus (remainder) function corner cases
Hi All, I'm curious as to what folks think about some choices in the compution of the remainder function. As an example where different choices can be made In [2]: -1e-64 % 1. Out[2]: 1.0 In [3]: float64(-1e-64) % 1. Out[3]: 0.99999999999999989 The first is Python, the second is in my branch. The first is more accurate numerically, but the modulus is of the same magnitude as the divisor. The second maintains the convention that the result must have smaller magnitude than the divisor. There are other corner cases along the same lines. So the question is, which is more desirable: maintaining numerical accuracy or enforcing mathematical convention? The differences are on the order of an ulp, but there will be a skew in the distribution of the errors if convention is maintained. 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 ;) Chuck
On Sat, Feb 13, 2016 at 9:31 AM, Charles R Harris <charlesr.harris@gmail.com
wrote:
Note that the other enforced bound is that the result have the same sign as the divisor. Python enforces that by adjusting the integer part, I enforce it by adjusting the remainder. Chuck
2016-02-13 17:42 GMT+01:00 Charles R Harris <charlesr.harris@gmail.com>:
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)).
On Sun, Feb 14, 2016 at 12:35 PM, Nils Becker <nilsc.becker@gmail.com> wrote:
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
? <snip> Chuck
On Sun, Feb 14, 2016 at 1:11 PM, Charles R Harris <charlesr.harris@gmail.com
wrote:
However, I do note that some languages offer two versions of modulus, one floor based and the other trunc based (effectively fmod). What I wanted is to keep the remainder consistent with the floor function in the C library. Chuck
On Sat, Feb 13, 2016 at 9:31 AM, Charles R Harris <charlesr.harris@gmail.com
wrote:
Note that the other enforced bound is that the result have the same sign as the divisor. Python enforces that by adjusting the integer part, I enforce it by adjusting the remainder. Chuck
2016-02-13 17:42 GMT+01:00 Charles R Harris <charlesr.harris@gmail.com>:
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)).
On Sun, Feb 14, 2016 at 12:35 PM, Nils Becker <nilsc.becker@gmail.com> wrote:
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
? <snip> Chuck
On Sun, Feb 14, 2016 at 1:11 PM, Charles R Harris <charlesr.harris@gmail.com
wrote:
However, I do note that some languages offer two versions of modulus, one floor based and the other trunc based (effectively fmod). What I wanted is to keep the remainder consistent with the floor function in the C library. Chuck
participants (2)
-
Charles R Harris -
Nils Becker