Overflow in linspace construction

[timholy]

Filing this so I don't forget to tackle it when I have time (not that I mind if someone beats me to it 😉):

julia> linspace(Float16(0.1094), Float16(0.9697), 300)
NaN16:Float16(0.002878):NaN16

Diagnosis

The problem is that these product+divisions compute the product first in TwicePrecision{Float16} before computing any divisions. This overflows the maximum allowed Float16.

One solution might be to interleave the products and divisions, but this appears to lead to ulp-level inaccuracies on other problems (I did not spend time chasing down why).

What those lines, plus the next, are actually doing is computing

((len-imin) * start_n + (imin-1) * stop_n) / (den * (len-1))

in a manner resistant to integer overflow (but fails to take into account floating point overflow).

Potential solution

Right now TwicePrecision implicitly is focused on AbstractFloat types, but perhaps we should also define methods appropriate for TwicePrecision{<:Integer}. If we implemented * and +, then we could compute that numerator without roundoff error. It seems theoretically possible that this might also contribute to fixing #20521.

This would be a step towards making TwicePrecision objects first-class, rather than their current usage which is solely intended for supporting Range computations. (Currently they are not even TwicePrecision<:Number, to avoid risk of ambiguities.) It's a bit questionable whether one wants to commit to broadly supporting arithmetic with TwicePrecision, though, so an alternative might be to implement the necessary splitprec, mul2, and add2 operations "manually" within that function body.

EDIT: or we could do that multiply/add in Int128 precision. That seems much easier. The conversion to type T is still delicate, however.