Skip to content
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

RFC: RoundNearest argument for rem #10946

Merged
merged 7 commits into from Jan 20, 2017
Merged

RFC: RoundNearest argument for rem #10946

merged 7 commits into from Jan 20, 2017

Conversation

simonbyrne
Copy link
Contributor

This is my attempt at resolving #9283. It adds an extra argument to rem, so that

rem(x,y,r::RoundingMode) = x - y*round(x/y,r)

without intermediate rounding, so that the result is always in the interval [-abs(y)/2,abs(y)/2]. This corresponds to the remainder function in the IEEE754 spec.

I've also implemented the corresponding rem2pi function, which is more accurate, and seems to be more useful, than mod2pi (I actually replaced the only 2 occurrences of mod2pi in Base, as this is what they were actually trying to get at).

If we're happy with this, I can add some docs, tests, and implement corresponding div functions.

@simonbyrne
Copy link
Contributor Author

And I guess this would also be useful to have for integer arguments.

@simonbyrne
Copy link
Contributor Author

Any suggestions for why the build is breaking on 32bit?

@pao
Copy link
Member

pao commented Apr 22, 2015

Tuple overhaul related. Should be fixed by 7ec501d. Either force-push or get someone who can rerun the tests to trigger it and you should be good to go.

@JeffBezanson
Copy link
Sponsor Member

Looks really good; please continue!

@StefanKarpinski
Copy link
Sponsor Member

Yes, this is great. I like the design.

@StefanKarpinski StefanKarpinski added the maths Mathematical functions label Sep 13, 2016
@StefanKarpinski StefanKarpinski added this to the 0.6.0 milestone Sep 13, 2016
@StefanKarpinski
Copy link
Sponsor Member

@simonbyrne – this was generally well-liked, can we get it rebased and merged?

@@ -397,6 +397,8 @@ export
reim,
reinterpret,
rem,
rem1,
rem2pi,
Copy link
Contributor

Choose a reason for hiding this comment

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

I thought we removed rem1, and does rem2pi need to be exported?

Copy link
Sponsor Member

Choose a reason for hiding this comment

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

🤷‍♂️

Copy link
Sponsor Member

Choose a reason for hiding this comment

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

Honestly, I'm not sure. I'm also not sure why rem2pi takes a rounding argument but mod2pi doesn't.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

mod2pi(x) is rem2pi(x, RoundDown)

@StefanKarpinski
Copy link
Sponsor Member

Ok, I resolved the conflicts as best I could, but to be perfectly honest, I have no idea what I'm doing here. Some are pretty clear (docstring added on master right below a new definition in this PR – keep both), others were way over my head – e.g. the previous modification to logdet in base/linalg/lu.jl. cc: @simonbyrne, @KristofferC, @andreasnoack, @jiahao, @sarvjeetsinghghotra.

airyai, airyaiprime, airybi, airybiprime,
airyaix, airyaiprimex, airybix, airybiprimex,
clamp, modf, ^, mod2pi, rem2pi,
airy, airyai, airyprime, airyaiprime, airybi, airybiprime, airyx,
Copy link
Contributor

Choose a reason for hiding this comment

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

wrong conflict resolution on the airy names I think

rem2pi(x::Float32, ::RoundingMode{:Nearest}) = Float32(rem2pi(Float64(x),RoundNearest))
rem2pi(x::Int32, ::RoundingMode{:Nearest}) = rem2pi(Float64(x),RoundNearest)
function rem2pi(x::Int64, ::RoundingMode{:Nearest})
fx = Float64(x)
Copy link
Contributor

Choose a reason for hiding this comment

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

4 space indent

@tkelman tkelman added needs docs Documentation for this change is required needs tests Unit tests are required for this change labels Jan 4, 2017
@simonbyrne
Copy link
Contributor Author

Updated, but still needs some tests.

`[-abs(y)/2, abs(y)/2]`.

- if `r == RoundToZero` (default), then the result is exact, and in the interval
`[0,abs(y))` if `x` is positive, or `(abs(y),0]` otherwise.
Copy link
Contributor

Choose a reason for hiding this comment

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

maybe math-mode (double backtick) the intervals since that isn't Julia notation?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Good call.

@simonbyrne
Copy link
Contributor Author

Two remaining questions:

  • should we deprecate mod2pi?
  • should we give rem2pi a default rounding mode?

x - 2π*round(x/(2π),r)

without any intermediate rounding. This internally uses a high precision approximation of
2π, and so will give a more `rem(x,2π,r)`
Copy link
Contributor

Choose a reason for hiding this comment

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

a more ... accurate result than ?

be the most accurate result.

- if `r == RoundToZero`, then the result is in the interval `[0,2π]` if `x` is positive,.
or `[2π,0]` otherwise.
Copy link
Contributor

Choose a reason for hiding this comment

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

[-2π,0] ?

@tkelman tkelman removed the needs docs Documentation for this change is required label Jan 5, 2017
@StefanKarpinski
Copy link
Sponsor Member

Why is rem2pi preferable to mod2pi? Doesn't one need each in different situations? Shouldn't they both take rounding arguments and default to the current rounding mode?

0.7853981633974481
```
"""
mod2pi(x) = rem2pi(x,RoundDown)
Copy link
Contributor Author

Choose a reason for hiding this comment

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

@StefanKarpinski mod2pi is now just a special case of rem2pi

@simonbyrne
Copy link
Contributor Author

simonbyrne commented Jan 9, 2017

I've added some tests. Comments welcome.

@tkelman tkelman removed the needs tests Unit tests are required for this change label Jan 9, 2017
@StefanKarpinski
Copy link
Sponsor Member

Can't we just make rem2pi a method of mod2pi then?

@simonbyrne
Copy link
Contributor Author

I think rem is more general that mod (mod to me implies that the values are always in the interval [0,x) for some x > 0).

@tkelman
Copy link
Contributor

tkelman commented Jan 19, 2017

rebase then merge IMO?

@simonbyrne simonbyrne merged commit cb1a47e into master Jan 20, 2017
@simonbyrne simonbyrne deleted the sb/rem-nearest branch January 20, 2017 19:52
@tkelman tkelman added the needs news A NEWS entry is required for this change label Jan 20, 2017
@tkelman
Copy link
Contributor

tkelman commented Jan 20, 2017

might be worth adding a news entry sometime before the release

@StefanKarpinski
Copy link
Sponsor Member

That is not correct: mod means that the value is in [0,m) for positive m or (m,0] for negative m, I.e. that the result sign always matches the modulus; rem means that the sign matches the argument rather than the modulus. For mathematical purposes, one generally wants mod not rem.

@simonbyrne
Copy link
Contributor Author

I took rem to mean remainder after division, as described in #9283. In that case mod is a special case (being when the division is rounded down).

@StefanKarpinski
Copy link
Sponsor Member

@simonbyrne and I just had the following chat conversation which cleared this issue up for me.


Stefan Karpinski [2:40 PM]

@simonbyrne: I’m a bit unclear on what the difference between my understanding of rem and yours is

Simon Byrne [2:44 PM]

I'm not sure I understand what you're proposing
are you saying we should just make rem(x, y, r::RoundingMode), mod(x, y, r::RoundingMode)?

Stefan Karpinski [2:59 PM]

I guess I don’t understand why rem is more fundamental than mod in your definition of it
@simonbyrne ^
(sorry, trying to get this happening as a synchronous conversation)

Simon Byrne [3:00 PM]

to me, mod(x,y) implies that the result must be in the interval [0,y) (if y>0)

Stefan Karpinski [3:01 PM]

the more general definition is that mod(x,y) is in [0,y) for positive y and (y,0] for negative y

Simon Byrne [3:01 PM]

sure

Stefan Karpinski [3:01 PM]

which is usually what you actually want for a mathematical modulus operation

Simon Byrne [3:01 PM]

yes
but then what do you call the operation r = f(x,y) which gives -abs(y)/2 <= x <= abs(y)/2?
I tried to explain my thinking here:
#9283
The idea being that rem and mod are both the same sort of thing
finding r such that x = q*y + r for some integer q
and the difference is how you compute q
(or, really, the direction in which you round q)
and so "remainder" seems like the most general term for this operation
(similarly, we could also make fld and cld special cases of div(x,y,::RoundingMode))

Stefan Karpinski [3:06 PM]

Why make rem the fundamental one and not mod?

Simon Byrne [3:06 PM]

because you're finding a remainder
not a modulus

Stefan Karpinski [3:06 PM]

Can you explain how rounding is the difference between rem and mod?

Simon Byrne [3:08 PM]

so current rem is defined as
rem(x,y) = x - y*trunc(x/y) (edited)
(without intermediate rounding)
mod(x,y) = x - y*floor(x/y) (edited)
so #10946 defined
rem(x,y,r::RoundingMode) = x - y*round(x/y, r)
i.e. the remainder after Euclidean division with rounding direction r

Stefan Karpinski [3:11 PM]

I see
So you can get the operation you were talking about via a different rounding mode

Simon Byrne [3:12 PM]

yes, that would just be
rem(x,y, RoundNearest) = x - y*round(x/y)

Stefan Karpinski [3:13 PM]

standard rem, mod correspond to round towards zero (trump), and round down (floor), respectively

Simon Byrne [3:13 PM]

yes
exactly

Stefan Karpinski [3:13 PM]

Ok, I'm starting to get this now

Simon Byrne [3:14 PM]

my plan is also to define corresponding methods for div as well (edited)
since that is useful (e.g. for implementing fixed point arithmetic)

Stefan Karpinski [3:14 PM]

So the question becomes what function to make the three-argument form a method of

Simon Byrne [3:15 PM]

yes

Stefan Karpinski [3:15 PM]

A new function entirely, rem or mod

Simon Byrne [3:15 PM]

yes
personally, remainder seems like the best name for it

Stefan Karpinski [3:16 PM]

Is the issue with mod2pi that the behavior it implements does not correspond to what the mod function does?

Simon Byrne [3:17 PM]

no
our current mod2pi(x) = mod(x, 2pi)
which is okay
but not that useful

Stefan Karpinski [3:17 PM]

But with better precision

Simon Byrne [3:17 PM]

it's better to do the argument reduction to [-pi, pi], rather than [0, 2pi]

Stefan Karpinski [3:18 PM]

Ah, I see

Simon Byrne [3:18 PM]

e.g. if you get x = -0.000001

Stefan Karpinski [3:19 PM]

The trouble is that this operation has no standard name

Simon Byrne [3:19 PM]

and so that corresponds to rem(x, 2pi, RoundNearest)

Stefan Karpinski [3:19 PM]

Right, right

Simon Byrne [3:19 PM]

hence rem2pi(x, RoundNearest)
i mean we could punt on dispatch

Stefan Karpinski [3:19 PM]

?

Simon Byrne [3:19 PM]

define an twopi::Irrational
and rem(x, twopi, RoundNearest)

Stefan Karpinski [3:20 PM]

I think you mean tau

Simon Byrne [3:20 PM]

sure

Stefan Karpinski [3:21 PM]

I did add tau in a PR at one point, but I think I'd generally rather not start dispatching on irrationals too much

Simon Byrne [3:21 PM]

yeah

Stefan Karpinski [3:21 PM]

Gives people the wrong idea

Simon Byrne [3:21 PM]

it's not that useful
and could trigger a lot of unnecessary bikeshedding
best to leave that to Tau.jl

Stefan Karpinski [3:22 PM]

If pi works, why doesn't 2pi, and 3pi, and pi/2, etc
Sure

Simon Byrne [3:22 PM]

yeah
C calls rem(x,y,RoundNearest) remainder (edited)

Stefan Karpinski [3:23 PM]

I could actually imagine extending irrational to having integer and rational multiples

Simon Byrne [3:23 PM]

as does the IEEE spec

Stefan Karpinski [3:23 PM]

But that's way too big a change to drag into this

Simon Byrne [3:23 PM]

yeah
agreed
maybe for 2.0
or 3.14 if we're going for TeX numbering

Stefan Karpinski [3:24 PM]

Is that how rem is defined for floating point in C?

Simon Byrne [3:24 PM]

no rem and remainder are different

Stefan Karpinski [3:24 PM]

I.e. In a way that is incompatible with integer rem?
Oh, I see. It calls that operation remainder

Simon Byrne [3:25 PM]

C has fmod (which is rem(x::Float64, y::Float64)
and remainder (which is rem(x::Float64, y::Float64, RoundNearest))

Stefan Karpinski [3:26 PM]

I wonder if we should run this by a few people with good judgement on this sort of thing, e.g. @JeffBezanson and @stevengj

Simon Byrne [3:27 PM]

and remquo, which gives the same answer as remainder, but also the last few bits of the quotient via a pointer
it doesn't define a mod function

Stefan Karpinski [3:27 PM]

I can see rem being a good place for this but something feels oddly asymmetrical about it

Simon Byrne [3:27 PM]

which is why we had all the hassles defining it back in the day
(the day being 2013 or so)

Stefan Karpinski [3:28 PM]

Ok, so your plan is to make rem/div the fundamental pair with three-arg methods

Simon Byrne [3:28 PM]

yes

Stefan Karpinski [3:29 PM]

And then conceptually make all the two-arg versions aliases for calling those with various rounding modes

Simon Byrne [3:29 PM]

yes

Stefan Karpinski [3:29 PM]

What about rem2pi(x)
Which rounding more does that imply?

Simon Byrne [3:30 PM]

the same as rem(x, 2pi), i.e. RoundToZero

Stefan Karpinski [3:30 PM]

Sign trunc, i.e. Sign matches x?

Simon Byrne [3:30 PM]

yes

Stefan Karpinski [3:30 PM]

And round to nearest doesn't get a name

Simon Byrne [3:30 PM]

i mean, we could make it the default?
but I think people would find that confusing

Stefan Karpinski [3:31 PM]

Yes, that seems pretty confusing

Simon Byrne [3:31 PM]

also would break a lot

Stefan Karpinski [3:31 PM]

Ok, thanks for explaining this to me
You've sold me on your grand vision

Sacha0 added a commit to Sacha0/julia that referenced this pull request May 11, 2017
@tkelman tkelman removed the needs news A NEWS entry is required for this change label May 16, 2017
tkelman pushed a commit that referenced this pull request May 16, 2017
* Add NEWS.md entry for three-argument rem and the new rem2pi (#10946).

* Clarify precise behaviour
tkelman pushed a commit that referenced this pull request May 17, 2017
* Add NEWS.md entry for three-argument rem and the new rem2pi (#10946).

* Clarify precise behaviour

(cherry picked from commit 8f9a948)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
maths Mathematical functions
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

5 participants