Julia: Euclidean division: `div`, `rem` and friends

[simonbyrne]

For any two real numbers x and y, Euclidean division is the problem of finding an integer q (the quotient) and number r (the remainder, on an interval of length abs(y)) such that

x = q*y + r

We currently offer 3 methods for computing q and r

| quotient function | equivalent exact expression | remainder function |
| --- | --- | --- |
| div(x,y) | trunc(x/y) | rem(x,y) |
| fld(x,y) | floor(x/y) | mod(x,y) |
| cld(x,y) | ceil(x/y) | mod(x,-y) |

Unfortunately, none of these actually correspond to the remainder function in the IEEE-754 standard, which is defined as the remainder after x/y has been rounded to the nearest integer, ties to even (this is also the remainder function in C).

My proposal is then to implement these by dispatching div, rem and divrem on RoundingMode. This is similar to what we do for round. e.g.:

div(x,y, ::RoundingMode{:ToZero}) = div(x, y)
rem(x,y, ::RoundingMode{:ToZero}) = rem(x, y)
div(x,y, ::RoundingMode{:Down}) = fld(x, y)
rem(x,y, ::RoundingMode{:Down}) = mod(x, y)
div(x,y, ::RoundingMode{:Up}) = cld(x, y)
rem(x,y, ::RoundingMode{:Up}) = mod(x, -y)

as well as appropriate definitions for RoundNearest, RoundNearestTiesAway and RoundNearestTiesUp. These would all have the properties that:

• div(x, y, r) is equal to round(x/y, r) (if computed without intermediate rounding)
• x == div(x, y, r)*y + rem(x, y, r) (if computed without intermediate rounding)

The first steps toward this were made in #10946 (defining rem(x,y, ::RoundingMode) for floating point arguments, but the remainder (no pun intended) is still to be done.

Another thing that would be useful to have, is something like C's remquo function, which returns the congruent modulo of the quotient, which is handy for massive range reductions (e.g. for sind).

[simonbyrne]

How would people feel about having div always return an Integer? I had a look through Base, and it seems that in most cases it is subsequently converted to an Int anyway.

[ScottPJones]

Always in Int? What if my arguments were BigInt? Or I used a new decimal floating point type? (the result would be an integer value, but you don't want to lose that it is a decimal floating point type).

[ScottPJones]

@simonbyrne I think you are confusing the concrete type Int with the abstract concept of the result being an integer value... I think it works correctly now.

[tkelman]

For div, there is the case of div(1.0e250, 1), would that have to throw an OverflowError to remain type-stable?

[simonbyrne]

@ScottPJones

Always in Int? What if my arguments were BigInt?

I should have expanded more. What I was thinking is that div(::Float64,::Float64) would return an Int, and div(::BigFloat,::BigFloat) would return a BigInt (similar to convert(Integer,x)). Though we could define div{T<:Integer}(::Type{T},x,y) to denote return type (in fact, this might be a good place to start, even if we don't make the full change).

Or I used a new decimal floating point type? (the result would be an integer value, but you don't want to lose that it is a decimal floating point type).

I guess this was the crux of my question: are there any use cases where it would be more useful to return a value of the same type as the arguments, as opposed to a subtype of Integer?

@tkelman I was actually thinking div should "wrap-around" (similar to remquo, but up to the full precision of the integer type), so that div(1e250,1) == 0. But we could also have a checked_div as well.

[simonbyrne]

Though we might have to throw errors in the case of Infs and NaNs.

[ScottPJones]

@simonbyrne I know I would not like div to return a different type than what I gave it... I think it might mess up later calculations, and how would you know what the mapping should be? Float32 goes to Int32, Float64 -> Int64, Float128 -> Int128, maybe, but it still seems a bit of a mess that you have to know in div what Integer type can hold any possible integral value of an arbitrary type.

[simonbyrne]

No, the return types would just be Int (for non-big types), or BigInt (for BigFloat/BigInt), i.e. the same as convert(Integer,x).

The opposite question is then: if we return a floating point value of the same type, what should we return when the answer is greater than maxintfloat? As a simple example,

julia> div(Int(8*maxintfloat(Float64)),5) # exact value
14411518807585587

Now, 14411518807585587 is not exactly representable by a Float64, so we would have to round it to either 14411518807585586 or 14411518807585588. Under default round-to-nearest, ties-to-even rules, this would be 14411518807585588, which is greater than the actual true value.

[ScottPJones]

@simonbyrne What is the problem there? you are doing a div on an Int, not on a Float64...

[kmsquire]

julia> Int64(div(8*maxintfloat(Float64), 5))
14411518807585588

Off by one.

[ScottPJones]

I still don't see how that has anything to do with div...
If you do maxintfloat(Float64), to get a Float64.
If you then multiply it by 8, you lose bits...
If you want a correct answer, _before_ you multiply by 8, promote it to a larger type...
(and Int is NOT larger for integers... whenever you write Int, think _Int32_!)

i.e.

julia> Int(div(BigFloat(maxintfloat(Float64))*8,5))
14411518807585587

and this is what would happen on a 32-bit machine:

ulia> div(Int(8*maxintfloat(Float64)),5)
ERROR: InexactError()
 in call at base.jl:36

[kmsquire]

Not relevant to the question, but thanks for pointing out the Int vs Int64 problem. Corrected inline above.

8*maxintfloat(Float64) is exactly representable as a Float64, but you are correct in that it has lost bits, of course.

Anyway, FWIW, I agree that div returning the same type as the operands is the least surprising behavior.

[ScottPJones]

It is only representable because you picked a number where it increased the exponent... but you are definitely past the range of representable integers.
I guess the problem here is that Julia normally doesn’t automatically promote from one type to a larger one, does it?

[kmsquire]

Well, Simon picked it, but yes, it was picked purposely. It's also perfectly representable (but the 15 numbers before or after are not).

And no Julia doesn't promote automatically, because that would preclude various useful optimizations.

Cheers!

[ScottPJones]

@kmsquire I wasn’t trying to say that Julia should promote automatically, but rather, that that is what provokes this “problem”. I still don’t see it as really being a problem though, I think if you care about exact results, you need to think about promoting things yourself, since Julia won’t do it for you.

[ScottPJones]

When I wrote a decimal floating point package ages ago, I simply “promoted” things internally to 128 bits before performing the operations (which were only 6, +, -, *, /, (div), and # (rem)). [the numbers were represented as 64-bit signed #s with a signed byte exponent]. The nice thing was that even on 16-bit platforms, the PDP-11, DG, and MS-DOS, you’d get exactly the same answer as on 32-bit machines, all out to 19 digits... (and this was back when there was no standard for binary floating point... you’d get different answers on every platform!). Could you not promote things just internally to the div() operation?

[simonbyrne]

Could you not promote things just internally to the div() operation?

What do you mean?

I'll leave div as it is for the moment, though I might define some div{T<:Integer}(::Type{T},x::Float64,y::Float64) methods, as these could be useful in a couple of places (such as rationalize).

One other slightly crazy idea is to define a UInt3 type, as then we can have a fast divrem(::Type{UInt3},x::Float64,y::Float64), via the C libm remquo function, which can be handy for range reduction (e.g. sind, cosd). Though it's probably not worth it until we have stack-allocated Refs.

[StefanKarpinski]

How would a three-bit unsigned integer type help?

[simonbyrne]

Technically you only need 2 bits: if you do div(x,90) then it will tell which quadrant of the circle you're on, and you can compute the appropriate sind/cosd of the remainder.

I'm not sure why remquo gives you 3, but that also happens to be what the x87 FPREM/FPREM1 instructions return for free when calculating remainders.

[simonbyrne]

I've updated the proposal (I've intentionally left off computing the trailing bits of the quotient, as that discussion seemed to have gotten sidetracked).

[StefanKarpinski]

I'm going to purge a bunch of irrelevant conversation from this thread, which can be seen here for posterity: https://gist.github.com/StefanKarpinski/353b2f22b18b0bb3a74638784f72ebc7.

[denisrosset]

Small remark: the operation described above is not Euclidean division, but rather truncated division.

The Euclidean remainder of rational numbers, floating point numbers is always 0, as division is always possible for nonzero divisors.

See https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/divmodnote-letter.pdf

For the Euclidean division of integers, the remainder should always be nonnegative -- I've just been bitten by that when implementing computation of Hermite forms (but in another language).

[simonbyrne]

@denisrosset thanks for the link, it's an interesting read.

I used the term "Euclidean" to refer to any division that computes an integer quotient. I wanted to avoid "integer division", as that might imply "division of integers". Unfortunately it seems that computer scientists just seem to call it "division" ignoring completely ignoring numbers other than integers. Avoiding confusing/ambiguous terminology was exactly the purpose of this issue.

Note that the Euclidean remainder can be non-zero for floating point numbers, since the quotient is required to be an integer, i.e.

julia> rem(1.5,1.0)
0.5

[denisrosset]

@simonbyrne : on an Euclidean domain, both q and r are members of the Euclidean domain.

• Ring of univariate polynomials over the reals: we have 3x^3 + 1 divided by 2x + 1 gives (3/2 x^2 - 3/4 x + 3/8)(2x + 1) + 5/8.

• Fields: all nonzero elements have an inverse, so q = x/y and r=0 always (for nonzero y).

What we have here is truncated division, which needs only an additive group and a total order, and indeed works on floating point numbers, rational numbers, algebraic numbers, etc... as you mention.

[simonbyrne]

Closed by #33040.

[StefanKarpinski]

Closing a four-digit issue ain't bad these days 🎉