Julia: revert linspace to creating a vector

[StefanKarpinski]

See https://github.com/JuliaLang/julia/issues/14420 and various other issues. It seems best to remove the LinSpace type from Base and move it to a package, keeping the linspace function but having it just return a vector of values.

[tkelman]

We already have several spellings for getting the linspace result as a vector, either collect or an equivalent comprehension. We wouldn't have a spelling for getting the lightweight iterator object of a range with exact start, end, and number of items if we did this. Tim's suggestion at https://github.com/JuliaLang/julia/pull/17408#issuecomment-234806131 for simplifying its behavior seemed popular so far.

[StefanKarpinski]

That was originally my thinking as well, but currently we don't have an easy spelling for getting linspace result as a vector with always correct end-points.

[StefanKarpinski]

I.e. https://github.com/JuliaLang/julia/issues/14420. @timholy's suggestion is certainly appealingly simple, but it really whiffs on a huge number of cases for internal points. I really don't think it makes sense to have this type and have it return half-baked results. Especially when the FloatRange behavior would actually be what one wants in those same cases.

[mbauman]

I disagree. If you care about the steps looking right, use the syntax that specifies a step. If you care about the number of elements and final endpoint, use the syntax that specifies those.

They're different algorithms. With different spellings. The only thing I'd change about Tim's implementation is how it prints itself; it should look different from a FloatRange, because it is.

Let's keep it as a lazy type. Folks are going to get more accustomed to all the specialized array types as we start using them more. And support for BLAS-incompatible array types is only getting better.

[StefanKarpinski]

export broken_linspace and correct_linspace? :P

[timholy]

I think my brain spends 10% of its cycles running @test_approx_eq_eps x y factor_of_two, so when I look at 0.3 and 0.299999999 I have to rub my eyes thrice to even notice that there's a difference.

Jokes aside, I do think there's a legitimate issue of whether it's really worth calling such differences "broken." As you yourself have pointed out, there are limits to representing mathematical quantities with a finite number of bits.

[eschnett]

IEEE goes to great lengths to specify that addition, multiplication, and a few other operations need to be calculated correctly up to the last bit. With @fastmath, I like to take shortcuts, but without, I expect all the precision I can get.

linspace is mentally a sufficiently simple operation that it should work without rounding error.

[StefanKarpinski]

... and for that you need do use the division operator.

[timholy]

Referencing https://github.com/JuliaLang/julia/pull/17408#issuecomment-234806131, this version seems likely to get it essentially perfect:

function Base.getindex(r::NewLinSpace, i::Integer)
    d = r.len-1
    s1 = fraction(r.len-i, r.leninv, d)
    s2 = fraction(i-1,     r.leninv, d)
    s1*r.start + s2*r.stop
end

@inline function fraction(k, dinv, d)
    x = k*dinv
    x + (k-x*d)*dinv
end

The only problem is it's 50% slower than our current LinSpace on a sum-over-all entries benchmark (eliding boundschecks). Perhaps there's something clever we could do to speed this up?

[StefanKarpinski]

Closed in favor of https://github.com/JuliaLang/julia/pull/18777.