Test of linear interpolation between floating point numbers

lerptest.jl

function lerp0(j, d, a, b)
    linspace(a,b,d+1)[j+1]
end
function lerp1(j, d, a, b)
    s = (d-j)/d; t = j/d
    a*s + b*t
end
function lerp2(j, d, a, b)
    (a*(d-j) + b*j)/d
end
function lerp3(j, d, a, b)
    # based on #18777, but make t < 0.5
    if j > div(d,2)
        j = d-j
        a,b = b,a
    end
    t = j/d
    fma(t, b, fma(-t, a, a))
end
function lerp4(j, d, a, b)
    # original #18777 proposal
    if j <= div(d,2)
        j = d-j
        a,b = b,a
    end
    t = j/d
    fma(t, b, fma(-t, a, a))
end


function lerp5(j, d, a, b)
    if a > b
        j = d-j
        a,b = b,a
    end
    a + (j/d)*(b-a)
end

@inline trunclo(x::Float64) = reinterpret(Float64, reinterpret(UInt64, x) & 0xffff_ffff_f800_0000)

function lerp8(j, d, a, b)
    if a > b
        j = d-j
        a,b = b,a
    end
    id_hi = trunclo(1/d)
    id_lo = fma(id_hi,-d,1)/d
    
    u = id_hi*j
    v = id_lo*j
    w = u + v
    
    wt = trunclo(w)
    wl = (u-wt)+v

    x = b-a
    xt = trunclo(x)
    xl = abs(a) > abs(b) ? b-(xt+a) : (b-xt)-a
    
    y = w*x
    yl = (((wt*xt-y)+wt*xl)+wl*xt)+wl*xl

    z = a+y
    zl = abs(a) > abs(y) ? (a-z)+y : (y-z)+a
    z + (yl+zl)
end


lerp_big(j,d,a,b) = lerp2(big(j),big(d),big(a),big(b))




function testlerp(f,n=1_000_000)
    u = 0.0
    Xu = (0,0,0.0,0.0)

    for i = 1:n
        a = randexp()
        b = randexp()
        j = rand(1:1000)
        k = rand(1:1000)
        d = j+k

        X = j,d,a,b

        x = f(X...)
        y = lerp_big(X...)
        r = oftype(x,abs(x-y)/eps(oftype(x,y)))
        if r > u
            Xu = X
            u = r
        end
    end
    u, Xu
end

@show testlerp(lerp0)
@show testlerp(lerp1)
@show testlerp(lerp2)
@show testlerp(lerp3)
@show testlerp(lerp4)
@show testlerp(lerp5)
@show testlerp(lerp8)

For lerp8: one problem with computing x = b-a is that it will run into trouble for the b = -a = realmax() linspace test. Discovered through hard experience 😉.