Julia: Range arithmetic not working with BigFloat

[ettersi]

julia> BigFloat(2)*(1:2)
ERROR: StackOverflowError:
Stacktrace:
 [1] range(::BigFloat, ::BigFloat, ::Int64) at ./range.jl:74 (repeats 80000 times)

julia> versioninfo()
Julia Version 0.6.2
Commit d386e40c17 (2017-12-13 18:08 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-2620M CPU @ 2.70GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, sandybridge)

[dlfivefifty]

Here's another one on Julia v0.7-beta:

julia> versioninfo()
Julia Version 0.7.0-beta.105
Commit 656d58782c* (2018-07-01 09:04 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin17.6.0)
  CPU: Intel(R) Core(TM) i7-7700HQ CPU @ 2.80GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, skylake)
Environment:
  JULIA_VERSION = 0.6

julia> (2:3) .- one(BigFloat)
ERROR: StackOverflowError:
Stacktrace:
 [1] _range(::BigFloat, ::BigFloat, ::Nothing, ::Int64) at ./range.jl:83 (repeats 80000 times)

[dlfivefifty]

This definition seems flawed:

# range.jl:86
_range(a::AbstractFloat, st::AbstractFloat, ::Nothing, len::Integer) = _range(promote(a, st)..., nothing, len)

It looks like there needs to be a

_range(a::T, st::T, ::Nothing, len::Integer) where T<:AbstractFloat = ?

I'm not sure what the default definition. I wonder if the definition for Float32, Float64 would work as is:

# twiceprecision.jl:424
function _range(a::T, st::T, ::Nothing, len::Integer) where T<:Union{Float16,Float32,Float64}
    start_n, start_d = rat(a)
    step_n, step_d = rat(st)
    if start_d != 0 && step_d != 0 &&
            T(start_n/start_d) == a && T(step_n/step_d) == st
        den = lcm(start_d, step_d)
        m = maxintfloat(T, Int)
        if abs(den*a) <= m && abs(den*st) <= m &&
                rem(den, start_d) == 0 && rem(den, step_d) == 0
            start_n = round(Int, den*a)
            step_n = round(Int, den*st)
            return floatrange(T, start_n, step_n, len, den)
        end
    end
    steprangelen_hp(T, a, st, 0, len, 1)
end

[dlfivefifty]

Even in 1.0 this is still broken:

julia> range(BigFloat(1), length=10)
ERROR: StackOverflowError:
Stacktrace:
 [1] _range(::BigFloat, ::BigFloat, ::Nothing, ::Int64) at ./range.jl:86 (repeats 80000 times)

[dlfivefifty]

@mbauman @andreasnoack Not sure who is in charge of ranges but it looks like this issue was missed. Can it get a 1.0.x target?

[andreasnoack]

I'm trying to avoid float range discussions.

[StefanKarpinski]

I can take a look at this, @timholy is the most recent float-range code author. Andreas disapproves of the whole thing so he's not going to get involved.

[timholy]

And @timholy got involved only through a desire to improve/extend the API but couldn't do that without diving into the whole high-precision range topic. Having spent far too much time on that, he confesses to being more than a little tired of the whole subject and would be grateful if others can dive in.

[dlfivefifty]

I don't actually care what it does, as long as it doesn't stack overflow. (Well, it would be nice if it gave something reasonably sensible.)

Could we just restore the previous implementation as the default for general AbstractFloat, and leave the more refined new implementation for Float64/Float32?

[timholy]

Yes, it's fine (aka, desirable) to bypass all the twiceprecision stuff for anything other than Float16/Float32/Float64.

[daanhb]

I'm afraid the issue lives on, sorry for intervening with #30944 and #30952, and thanks @dlfivefifty for pointing it out. In the latter, admittedly fairly trivial, pull request I had automatically assumed that sidestepping the twiceprecision technique was the right (safest) default thing to do.

(For the record, I quite like the twiceprecision technique and was planning on replacing my own crude implementation of equispaced grid points elsewhere (in an unregistered package) with it.)

[felixcremer]

This is closed by https://github.com/JuliaLang/julia/pull/30952.