solve_linspace_divisor.jl

solve_divisor{T}(::Type{T},a,b,N) = solve_divisor(convert(T, a), convert(T, b), convert(T, N))

function solve_divisor{T}(a::T, b::T, N::T)
    #x = nextfloat(typemin(T))
    x = T(1)
    while x < 2 #prevfloat(typemax(T))
        if ((a*x)*N)/(x*N) == a && ((b*x)*N)/(x*N) == b
            return x
        end
        x = nextfloat(x)
    end
    error("No workable value of x found; a = $a, b = $b, N = $N")
end

# Test all floating point intervals (a,b) deterministically by iterating
# through the sequence
function testall{T}(::Type{T}, N)
    a = T(0.12345679123456) #nextfloat(typemin(T))
    while true
        @show a
        b = T(1.12345679123456) #nextfloat(typemin(T))
        while true
            # Ignore cases which will overflow the range
            if abs(2*Float64(N)*a) < prevfloat(typemax(T)) &&
               abs(2*Float64(N)*b) < prevfloat(typemax(T))
                x = solve_divisor(T, a,b,N)
            end
            b = nextfloat(b)
            if b >= 10
                break
            end
        end
        a = nextfloat(a)
        if a >= 10
            break
        end
    end
end

function testall_random{T}(::Type{T})
    iter = 0
    while true
        solve_divisor(Float32, 10*rand(Float32), 10*rand(Float32), rand(1:100))
        iter += 1
        if iter % 1_000_000 == 0
            info("Iteration $iter")
        end
    end
end

# Exhaustive tests for Float16
#testall(Float16, 49)
#testall(Float16, 3)

# Following appears to be a counter example in Float16, needs more checking
# solve_divisor(Float16, 0.00012195, 1.999, 3)

#testall_random(Float32)

testall(Float32, 3)