More accurate Base.lerpi for higher-precision numbers

[ffevotte]

This should fix #37276.

Examples of things that are more accurate with this PR include:

julia> LinRange(Complex{BigFloat}(0), Complex{BigFloat}(1), 11)[2]
0.1000000000000000000000000000000000000000000000000000000000000000000000000000002 + 0.0im

julia> LinRange([0//1], [1//1], 11)[2]
1-element Vector{Rational{Int64}}:
 1//10

in comparison to the current master:

julia> LinRange(Complex{BigFloat}(0), Complex{BigFloat}(1), 11)[2]
0.1000000000000000055511151231257827021181583404541015625 + 0.0im

julia> LinRange([0//1], [1//1], 11)[2]
1-element Array{Rational{Int64},1}:
 3602879701896397//36028797018963968

I think that the proposed implementation of a generic Base.lerpi makes it less useful to have a specific method for rationals (such as defined in rational.jl:467).

One potential downside to the way things are done in this PR, is that there could potentially be a loss of accuracy w.r.t how things are currently done, when working on ranges of small (i.e. less than 64 bits) floating-point numbers. I'm not sure whether that was intended, but in the current implementation, the j/d quotient usually results in t being a Float64, meaning that computations involving a and b will be promoted to Float64 too if the end points are smaller FP numbers. With the current implementation, all calculations use the same precision as the range end points.

It would not be too difficult to restore the current master behavior for smaller FP number types, but I'm not sure whether it is worth making the implementation more complicated. In particular, I have no idea how many users expect and rely on the fact that points in a LinRange{Float16} are computed in double precision before being rounded to half precision. Would somebody have some insight on this?

[simeonschaub]

Is this perhaps a bit too loose? Would it make sense to set a lower tolerance here?

[ffevotte]

Thanks for the review!

Are you referring to the tolerance of isapprox?
AFAIU, the default tolerance used by isapprox is given by the square root of eps of the types of the operands, which we know to be BigFloats here. In this case, this means a relative tolerance of approximately 4e-39 (which makes everything fail if any part of the computation is performed in double precision).

But you're right that it is much higher than the few ulps that we could expect; we could very well set this to something like 1e-75 and still be safe. TBH, the reason why I kept the default tolerance parameters is because this seems to be how it's done for all other approximate equality tests in test/ranges.jl.

TL;DR: I don't have any strong opinion about this. I can very well change this tolerance to any value deemed more appropriate.

[simeonschaub]

Yes, that's what I meant. I think it's better to be more conservative with tests like these, because it's always better to discover regressions early. Perhaps the other tests should be more strict as well, but that's probably a separate discussion.

[ffevotte]

OK. I agree that more strictness is always a good thing in tests. I amended to commit to include a relative tolerance of 10eps(BigFloat), which I guess makes the intent explicit (although in this particular case results are exactly rounded, I think it is more future-proof to not rely on this and allow an error of a few ulps, which should correspond to the expectation users should have in the general case)

[JeffBezanson]

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

[ffevotte]

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

[timholy]

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

[timholy]

On balance perhaps the rational approach makes the most sense.

[ffevotte]

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

[timholy]

@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

[timholy]

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

[ffevotte]

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.

[timholy]

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

[ffevotte]

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.

[StefanKarpinski]

I believe @timholy originally introduced this for accurate linear interpolation of ranges, so he may have some idea.

[ffevotte]

I believe @timholy originally introduced this for accurate linear interpolation of ranges, so he may have some idea.

Thanks. I indeed found #18777 where Base.lerpi seems to have appeared. The entire discussion is well worth a read and very relevant here but, although a few (sophisticated) implementations of lerpi are discussed there, I haven't found much about the generic implementation that was chosen in the end.

I'd love to hear about the details if someone happens to remember them.

[timholy]

In particular, I have no idea how many users expect and rely on the fact that points in a LinRange{Float16} are computed in double precision before being rounded to half precision. Would somebody have some insight on this?

I don't have any useful insight. I doubt it would be a problem but I am not sure.