10000 BUG: Potential fix for divmod(1.0, 0.0) to raise divbyzero and return inf by anirudh2290 · Pull Request #16161 · numpy/numpy · GitHub
[go: up one dir, main page]

Skip to content

BUG: Potential fix for divmod(1.0, 0.0) to raise divbyzero and return inf #16161

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 6 commits into from
Nov 25, 2020
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Next Next commit
BUG: Potential fix for 1//0 to raise divbyzero and return inf
  • Loading branch information
anirudh2290 committed Jun 5, 2020
commit 64e525d4333e5cac8d04debfe205cb779f0e4d68
40 changes: 39 additions & 1 deletion numpy/core/src/npymath/npy_math_internal.h.src
Original file line number Diff line number Diff line change
Expand Up @@ -624,6 +624,39 @@ NPY_INPLACE @type@ npy_logaddexp2@c@(@type@ x, @type@ y)
}
}

/*
* Wrapper function for remainder edge cases
* Internally calls npy_divmod*
*/
NPY_INPLACE @type@
npy_remainder@c@(@type@ a, @type@ b)
{
@type@ mod;
if (!b) {
mod = npy_fmod@c@(a, b);
npy_set_floatstatus_invalid();
} else {
npy_divmod@c@(a, b, &mod);
}
return mod;
}

NPY_INPLACE @type@
npy_floor_divide@c@(@type@ a, @type@ b) {
@type@ div, mod;
if (!b) {
div = a / b;
if (!a) {
npy_set_floatstatus_invalid();
} else {
npy_set_floatstatus_divbyzero();
}
} else {
div = npy_divmod@c@(a, b, &mod);
}
return div;
}

/*
* Python version of divmod.
*
Expand All @@ -637,9 +670,14 @@ npy_divmod@c@(@type@ a, @type@ b, @type@ *modulus)
mod = npy_fmod@c@(a, b);

if (!b) {
div = a / b;
npy_set_floatstatus_invalid();
if (a) {
npy_set_floatstatus_divbyzero();
}
/* If b == 0, return result of fmod. For IEEE is nan */
*modulus = mod;
return mod;
return div;
}

/* a - mod should be very nearly an integer multiple of b */
Expand Down
20 changes: 14 additions & 6 deletions numpy/core/src/umath/loops.c.src
Original file line number Diff line number Diff line change
Expand Up @@ -2004,8 +2004,7 @@ NPY_NO_EXPORT void
BINARY_LOOP {
const @type@ in1 = *(@type@ *)ip1;
const @type@ in2 = *(@type@ *)ip2;
@type@ mod;
*((@type@ *)op1) = npy_divmod@c@(in1, in2, &mod);
*((@type@ *)op1) = npy_floor_divide@c@(in1, in2);
}
}

Expand All @@ -2015,7 +2014,7 @@ NPY_NO_EXPORT void
BINARY_LOOP {
const @type@ in1 = *(@type@ *)ip1;
const @type@ in2 = *(@type@ *)ip2;
npy_divmod@c@(in1, in2, (@type@ *)op1);
*((@type@ *) op1) = npy_remainder@c@(in1, in2);
}
}

Expand Down Expand Up @@ -2360,8 +2359,13 @@ HALF_floor_divide(char **args, npy_intp const *dimensions, npy_intp const *steps
BINARY_LOOP {
const npy_half in1 = *(npy_half *)ip1;
const npy_half in2 = *(npy_half *)ip2;
npy_half mod;
*((npy_half *)op1) = npy_half_divmod(in1, in2, &mod);

float fh1 = npy_half_to_float(in1);
float fh2 = npy_half_to_float(in2);
float div;

div = npy_floor_dividef(fh1, fh2);
*((npy_half *)op1) = npy_float_to_half(div);
}
}

Expand All @@ -2371,7 +2375,11 @@ HALF_remainder(char **args, npy_intp const *dimensions, npy_intp const *steps, v
BINARY_LOOP {
const npy_half in1 = *(npy_half *)ip1;
const npy_half in2 = *(npy_half *)ip2;
npy_half_divmod(in1, in2, (npy_half *)op1);
float fh1 = npy_half_to_float(in1);
float fh2 = npy_half_to_float(in2);
float mod;
mod = npy_remainderf(fh1, fh2);
*((npy_half *)op1) = npy_float_to_half(mod);
}
}

Expand Down
12 changes: 10 additions & 2 deletions numpy/core/src/umath/scalarmath.c.src
Original file line number Diff line number Diff line change
Expand Up @@ -285,7 +285,11 @@ static void
@name@_ctype_floor_divide(@type@ a, @type@ b, @type@ *out) {
@type@ mod;

*out = npy_divmod@c@(a, b, &mod);
if (!b) {
*out = a / b;
} else {
*out = npy_divmod@c@(a, b, &mod);
}
}


Expand Down Expand Up @@ -318,7 +322,11 @@ static void
half_ctype_floor_divide(npy_half a, npy_half b, npy_half *out) {
npy_half mod;

*out = npy_half_divmod(a, b, &mod);
if (!b) {
*out = a / b;
} else {
*out = npy_half_divmod(a, b, &mod);
}
}


Expand Down
7 changes: 7 additions & 0 deletions numpy/core/tests/test_scalarmath.py
Original file line number Diff line number Diff line change
Expand Up @@ -276,6 +276,10 @@ def test_float_modulus_corner_cases(self):
# Check nans, inf
with suppress_warnings() as sup:
sup.filter(RuntimeWarning, "invalid value encountered in remainder")
sup.filter(RuntimeWarning, "divide by zero encountered in remainder")
sup.filter(RuntimeWarning, "divide by zero encountered in floor_divide")
sup.filter(RuntimeWarning, "divide by zero encountered in divmod")
sup.filter(RuntimeWarning, "invalid value encountered in divmod")
for dt in np.typecodes['Float']:
fone = np.array(1.0, dtype=dt)
fzer = np.array(0.0, dtype=dt)
Expand All @@ -290,6 +294,9 @@ def test_float_modulus_corner_cases(self):
assert_(np.isnan(rem), 'dt: %s' % dt)
rem = operator.mod(finf, fone)
assert_(np.isnan(rem), 'dt: %s' % dt)
for op in [floordiv_and_mod, divmod]:
div, mod = op(fone, fzer)
assert_(np.isinf(div)) and assert_(np.isnan(mod))

def test_inplace_floordiv_handling(self):
# issue gh-12927
Expand Down
52 changes: 51 additions & 1 deletion numpy/core/tests/test_umath.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
assert_, assert_equal, assert_raises, assert_raises_regex,
assert_array_equal, assert_almost_equal, assert_array_almost_equal,
assert_array_max_ulp, assert_allclose, assert_no_warnings, suppress_warnings,
_gen_alignment_data, assert_array_almost_equal_nulp
_gen_alignment_data, assert_array_almost_equal_nulp, assert_warns
)

def on_powerpc():
A3E2 Expand Down Expand Up @@ -293,6 +293,19 @@ def test_floor_division_signed_zero(self):
assert_equal(np.signbit(x//1), 0)
assert_equal(np.signbit((-x)//1), 1)

@pytest.mark.parametrize('dtype', np.typecodes['Float'])
def test_floor_division_corner_cases(self, dtype):
# test corner cases like 1.0//0.0 for errors and return vals
x = np.zeros(10, dtype=dtype)
y = np.ones(10, dtype=dtype)
# divide by zero error check
with np.errstate(divide='raise', invalid='ignore'):
assert_raises(FloatingPointError, np.floor_divide, y, x)
# verify 1.0//0.0 computations return inf
with np.errstate(divide='ignore'):
z = np.floor_divide(y, x)
assert_(np.isinf(z).all())

def floor_divide_and_remainder(x, y):
return (np.floor_divide(x, y), np.remainder(x, y))

Expand Down Expand Up @@ -366,6 +379,32 @@ def test_float_remainder_roundoff(self):
else:
assert_(b > rem >= 0, msg)

@pytest.mark.parametrize('dtype', np.typecodes['Float'])
def test_float_remainder_errors(self, dtype):
# Check valid errors raised for divmod and remainder
fzero = np.array(0.0, dtype=dtype)
fone = np.array(1.0, dtype=dtype)
finf = np.array(np.inf, dtype=dtype)
with np.errstate(invalid='raise'):
assert_raises(FloatingPointError, np.remainder, fone, fzero)
# since divmod is combination of both remainder and divide
# ops it will set both dividebyzero and invalid flags
with np.errstate(divide='raise', invalid='ignore'):
assert_raises(FloatingPointError, np.divmod, fone, fzero)
with np.errstate(divide='ignore', invalid='raise'):
assert_raises(FloatingPointError, np.divmod, fone, fzero)
with np.errstate(invalid='raise'):
assert_raises(FloatingPointError, np.divmod, fzero, fzero)
with np.errstate(invalid='raise'):
assert_raises(FloatingPointError, np.divmod, finf, finf)

def test_float_remainder_overflow(self):
a = np.finfo(np.float64).tiny
with np.errstate(over='ignore', invalid='ignore'):
div, mod = np.divmod(4, a)
np.isinf(div)
assert_(mod == 0)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we test that these warnings are given, or did that fail on some hardware?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

we ensure the errors invalid and dividebyzero are raised for divmod in line 370 to 373. Is this what you are looking for ?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, I guess I meant the specific case for tiny.


def test_float_remainder_corner_cases(self):
# Check remainder magnitude.
for dt in np.typecodes['Float']:
Expand All @@ -379,20 +418,31 @@ def test_float_remainder_corner_cases(self):
# Check nans, inf
with suppress_warnings() as sup:
sup.filter(RuntimeWarning, "invalid value encountered in remainder")
sup.filter(RuntimeWarning, "invalid value encountered in divmod")
sup.filter(RuntimeWarning, "divide by zero encountered in divmod")
for dt in np.typecodes['Float']:
fone = np.array(1.0, dtype=dt)
fzer = np.array(0.0, dtype=dt)
finf = np.array(np.inf, dtype=dt)
fnan = np.array(np.nan, dtype=dt)
rem = np.remainder(fone, fzer)
assert_(np.isnan(rem), 'dt: %s, rem: %s' % (dt, rem))
div, rem = np.divmod(fone, fzer)
assert_(np.isinf(div), 'dt: %s, div: %s' % (dt, rem))
assert_(np.isnan(rem), 'dt: %s, rem: %s' % (dt, rem))
# MSVC 2008 returns NaN here, so disable the check.
#rem = np.remainder(fone, finf)
#assert_(rem == fone, 'dt: %s, rem: %s' % (dt, rem))
rem = np.remainder(fone, fnan)
assert_(np.isnan(rem), 'dt: %s, rem: %s' % (dt, rem))
rem = np.remainder(finf, fone)
assert_(np.isnan(rem), 'dt: %s, rem: %s' % (dt, rem))
div, rem = np.divmod(fzer, fzer)
assert_(np.isnan(rem), 'dt: %s, rem: %s' % (dt, rem))
assert_(np.isnan(div), 'dt: %s, rem: %s' % (dt, rem))
div, rem = np.divmod(finf, finf)
assert_(np.isnan(div), 'dt: %s, rem: %s' % (dt, rem))
assert_(np.isnan(rem), 'dt: %s, rem: %s' % (dt, rem))


class TestCbrt:
Expand Down
0