Skip to content

More accurate Base.lerpi for higher-precision numbers #37281

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.

Sign up for GitHub

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

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

More accurate Base.lerpi for higher-precision numbers #37281

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
2 changes: 1 addition & 1 deletion base/range.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -685,7 +685,7 @@ end

function lerpi(j::Integer, d::Integer, a::T, b::T) where T
@_inline_meta
t = j/d
t = eltype(T)(j)/d
Copy link
Member

Choose a reason for hiding this comment

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

Why does this use eltype instead of just T(j)?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

This was required to satisfy the tests at https://github.com/JuliaLang/julia/blob/master/test/ranges.jl#L1283-L1290
In short, it ensures that things like this work:

julia> LinRange([0., 10.], [1., 12.], 4)
4-element LinRange{Array{Float64,1}}:
 [0.0, 10.0],[0.333333, 10.666667],[0.666667, 11.33333],[1.0, 12.0]

(I actually discovered such uses were possible when tests failed with an earlier version of my fix!)

Copy link
Member

@timholy timholy Aug 30, 2020 •
edited
Loading

Choose a reason for hiding this comment

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

A problem is that the current implementation works even for objects that don't define eltype(T): basically it just needs an algebra supporting addition and multiplication by reals.

One option would be t = j // d but I will wager you that there will be overflow test failures somewhere. Another option would be t = oftype(one(T), j)/d but amazingly even one(::Type{Vector{T}}) where T = one(T) isn't defined. (Why not? That seems crazy.)

Copy link
Member

Choose a reason for hiding this comment

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

On balance perhaps the rational approach makes the most sense.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

A problem is that the current implementation works even for objects that don't define eltype(T): basically it just needs an algebra supporting addition and multiplication by reals.

That's true. Are there any examples of such types in the stdlib? If we want to support such uses, I'd like to add a specific test case for them (but I'd like to avoid defining a new type only for this purpose if possible).

Copy link
Member

@timholy timholy Aug 31, 2020 •
edited
Loading

Choose a reason for hiding this comment

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

@StefanKarpinski has this annoying tendency to write tests that include the extreme ranges (joking, if he doesn't do that you can be sure that some pranksterconscientious developer will report it as an issue).

julia> start, stop = 1, typemax(Int)
(1, 9223372036854775807)

julia> t = 3//4
3//4

julia> (1 - t) * start + t * stop
ERROR: OverflowError: 3 * 9223372036854775807 overflowed for type Int64
Stacktrace:
 [1] throw_overflowerr_binaryop(op::Symbol, x::Int64, y::Int64)
   @ Base.Checked ./checked.jl:154
 [2] checked_mul
   @ ./checked.jl:288 [inlined]
 [3] *(x::Rational{Int64}, y::Int64)
   @ Base ./rational.jl:314
 [4] top-level scope
   @ REPL[3]:1

Copy link
Member

Choose a reason for hiding this comment

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

Interestingly, though, if we reversed the order of operation (divide first, multiply second) then it would be fine---though then you should check underflow for ranges like one going from eps(0.0) to 3*eps(0.0). If you can find a robust solution, it almost seems worth creating a SmallRatio type for this (example: https://github.com/timholy/Ratios.jl).

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Thanks, I was so focused on FP issues that hadn't thought about integer ranges.

One thing that we're able to handle in the current master is the following (which hits overflow issues with the j//d implementation as you mentioned):

julia> LinRange{Int}(1, 1+2^62, 5)
5-element LinRange{Int64}:
 1,1152921504606846976,2305843009213693952,3458764513820540928,4611686018427387904

(I'll note that it's currently not possible to handle LinRange{Int}(1, typemax(Int), 3) because lerpi uses FP numbers and typemax(Int) is not representable as a Float64)


In any case, in light of this discussion, I'll try and think a bit more before proposing a new implementation that hopefully meets all envisioned use cases.

Copy link
Member

Choose a reason for hiding this comment

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

One last issue to pay attention to is performance. I think it will be slower the more divisions you do.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes, that's what I meant to say in my message above:

There might be a small performance impact though, since lerpi would now perform two divisions: one to compute t, and one to compute 1-t. Is that a problem? Would you like me to benchmark this?

I'll try and find an adequate compromise between all concerns involved in this.

T((1-t)*a + t*b)
end

Expand Down
9 changes: 9 additions & 0 deletions test/ranges.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1326,6 +1326,15 @@ end
@test r[4] === 3.0+im
end

@testset "issue #37276" begin
r = range(Complex{BigFloat}(0), stop=Complex{BigFloat}(1), length=4)
@test isa(@inferred(r[2]), Complex{BigFloat})
@test r[1] == 0.
@test r[2] ≈ big"1"/3 rtol=10eps(BigFloat)
@test r[3] ≈ big"2"/3 rtol=10eps(BigFloat)
@test r[4] == 1.
end

# ambiguity between (:) methods (#20988)
struct NotReal; val; end
Base.:+(x, y::NotReal) = x + y.val
Expand Down