[BUG] np.mod with output outside the divisor range

[juntyr]

np.mod(np.cosh(QuadPrecision("-11357.216553474703894801348310092223")), QuadPrecision("-1.7976931345860068e+308"))

# prints QuadPrecision('1.145669756344548161919429759370211e+4898', backend='sleef')

np.mod should return a result with the same sign as the divisor, her negative, and it should be in range of the absolute value of the divisor.

The issue likely comes from rounding errors in

numpy-quaddtype/src/include/ops.hpp

Lines 645 to 648 in 156edad

	Sleef_quad quotient = Sleef_divq1_u05(*a, *b);
	Sleef_quad floored = Sleef_floorq1(quotient);
	Sleef_quad product = Sleef_mulq1_u05(floored, *b);
	Sleef_quad result = Sleef_subq1_u05(*a, product);

[SwayamInSync]

Thanks for catching this @juntyr I can confirm this bug and also it gets resolved via using fmod path.
I will update the code to use it and also later inspect which sub-operation in the raw formula magnifies the ULP.

[juntyr]

@SwayamInSync Here's my workaround using fmod (which has lower ULP since we don't combine div + mul + sub all with rounding errors):

out = np.fmod(p, q)
# correct the quadrant of the remainder
np.add(out, q, out=out, where=((p < 0) & (q > 0)))
np.add(out, q, out=out, where=((p > 0) & (q < 0)))
# ensure that |out| < |q|
np.fmod(out, q, out=out)
# ensure correct handling for finite p and infinite q
np.copyto(out, p, where=(np.isfinite(p) & np.isinf(q) & (np.sign(p) == np.sign(q))))
np.copyto(
    out,
    q,
    where=(np.isfinite(p) & np.isinf(q) & (np.sign(p) != np.sign(q)) & (p != 0)),
)
# ensure correct sign for NaNs and zeros
np.copysign(out, q, out=out, where=(np.isnan(out) | (out == 0)))
return out

I find https://en.wikipedia.org/wiki/Modulo#Variants_of_the_definition very helpful

[juntyr]

Btw @ngoldbaum what would the process be for numpy to expose modulo based on round-to-even division (remainder in C) as a ufunc and function?

[ngoldbaum]

I think you'd probably have an uphill battle getting it into NumPy unless there's a standard way to spell that operation.

[juntyr]

I think you'd probably have an uphill battle getting it into NumPy unless there's a standard way to spell that operation.

Aside from bikeshedding about the name, would there be other concerns? It’s a very useful modulo mode to have, it’s available through math.h (and SLEEF supports it as well), so the implementation should not be too difficult?

In terms of naming, it could be rint_remainder (or more verbosely round_ties_even_remainder, the current remainder/mod would be floor_remainder, similar to floor_divide, fmod would be trunc_remainder; the extra aliases could be added or not)

[juntyr]

Or follow Rust with rem_euclid -> rem_floor, rem_trunc, rem_rint, …

[ngoldbaum]

Well the issue isn't really utility, it's more that adding things to the NumPy API has ecosystem-wide ripple effects and we're generally pretty conservative about adding new API surface.

You'd need to follow the process, which starts with digging through the numpy mailing list and issue tracker to look for prior discussions about this topic.