9283.txt

Subject: [julia] Euclidean division: `div`, `rem` and friends (#9283)
------------------------

From: Simon Byrne <notifications@github.com>
Date: Tue, Dec 9, 2014 at 6:35 AM
To: JuliaLang/julia <julia@noreply.github.com>


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).

One idea would be implement these by dispatching div, rem and divrem on RoundingMode (similar to what I proposed in #8750 (comment).

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).

—
Reply to this email directly or view it on GitHub.


----------
From: Jeff Bezanson <notifications@github.com>
Date: Tue, Feb 3, 2015 at 8:51 AM
To: JuliaLang/julia <julia@noreply.github.com>


Could we find a better name for the "maths" label? Everything is math. Or we could delete the label.


----------
From: Simon Byrne <notifications@github.com>
Date: Sun, Apr 26, 2015 at 2:40 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


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.


----------
From: Scott P. Jones <notifications@github.com>
Date: Sun, Apr 26, 2015 at 4:16 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


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).


----------
From: Scott P. Jones <notifications@github.com>
Date: Sun, Apr 26, 2015 at 4:20 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


@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.


----------
From: Scott P. Jones <notifications@github.com>
Date: Sun, Apr 26, 2015 at 4:24 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


@JeffBezanson "Everything is math." Soooo not true!  :grinning: I might counter with "Everything is strings!" (having spent way too long maintaining a language where that was true).


----------
From: Tony Kelman <notifications@github.com>
Date: Sun, Apr 26, 2015 at 4:29 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


"Maths" is rather unspecific as an issue-tracking label, but I don't think it matters that much.

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


----------
From: Simon Byrne <notifications@github.com>
Date: Sun, Apr 26, 2015 at 5:00 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


@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.


----------
From: Simon Byrne <notifications@github.com>
Date: Sun, Apr 26, 2015 at 5:06 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


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


----------
From: Jeff Bezanson <notifications@github.com>
Date: Sun, Apr 26, 2015 at 7:58 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


So strings are outside the realm of mathematics? I daresay you have to go much further afield to find a candidate extra-mathematical topic, perhaps consciousness or morality, but certainly not strings.


----------
From: Scott P. Jones <notifications@github.com>
Date: Sun, Apr 26, 2015 at 9:05 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


Well, I have some physicist friends who'd tell me everything in the universe is strings :grinning:
I didn't say that strings were outside the realm of mathematics either.
However, in general computing, you can generally represent all of your numbers as strings, but not vice-versa... that's what I was getting at!


----------
From: Scott P. Jones <notifications@github.com>
Date: Sun, Apr 26, 2015 at 9:25 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


@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.


----------
From: Simon Byrne <notifications@github.com>
Date: Sun, Apr 26, 2015 at 10:58 PM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


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.


----------
From: Scott P. Jones <notifications@github.com>
Date: Mon, Apr 27, 2015 at 9:05 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


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


----------
From: Kevin Squire <notifications@github.com>
Date: Mon, Apr 27, 2015 at 10:12 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


julia> Int(div(8*maxintfloat(Float64), 5))
14411518807585588
Off by one.


----------
From: Scott P. Jones <notifications@github.com>
Date: Mon, Apr 27, 2015 at 11:34 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


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

----------
From: Kevin Squire <notifications@github.com>
Date: Mon, Apr 27, 2015 at 12:27 PM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


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.


----------
From: Scott P. Jones <notifications@github.com>
Date: Mon, Apr 27, 2015 at 5:15 PM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


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?


----------
From: Kevin Squire <notifications@github.com>
Date: Mon, Apr 27, 2015 at 6:51 PM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


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!

On Monday, April 27, 2015, Scott P. Jones <notifications@github.com> wrote:

> 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?
>
> —
> Reply to this email directly or view it on GitHub
> <https://github.com/JuliaLang/julia/issues/9283#issuecomment-96894913>.

----------
From: Scott P. Jones <notifications@github.com>
Date: Tue, Apr 28, 2015 at 12:36 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


@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.


----------
From: Scott P. Jones <notifications@github.com>
Date: Tue, Apr 28, 2015 at 12:43 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


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?


----------
From: Simon Byrne <notifications@github.com>
Date: Tue, Apr 28, 2015 at 1:22 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


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(::Type{Int},x,y) 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,y), 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.


----------
From: Simon Byrne <notifications@github.com>
Date: Tue, Apr 28, 2015 at 1:34 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>


Technically you only need 2 bits: if you do div(x,90) then it will tell which quadrant of the circle you're on, then 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.


----------
From: Simon Byrne <notifications@github.com>
Date: Mon, Jan 23, 2017 at 10:49 AM
To: JuliaLang/julia <julia@noreply.github.com>
Cc: Stefan Karpinski <stefan@karpinski.org>, Comment <comment@noreply.github.com>


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).

—
You are receiving this because you commented.
Reply to this email directly, view it on GitHub, or mute the thread.