You do get the numpy answer if you do the division last. With
xp = [0.0,3600.0]
yp = [25.7,26]
x = 3240
f = x/(xp[2]-xp[1])
instead of
julia> (1-f)*yp[1] + f*yp[2]
25.970000000000002
do
julia> ((xp[2] - x)*yp[1] + (x - xp[1])*yp[2])/(xp[2] - xp[1])
25.97
If all of these were integers, though, you’d be vulnerable to overflow:
xp = [0, 3600]
yp = [typemax(Int), typemax(Int)]
julia> ((xp[2] - x)*yp[1] + (x - xp[1])*yp[2])/(xp[2] - xp[1])
-1.0
but the f version is not:
julia> f = x/(xp[2]-xp[1])
0.9
julia> (1-f)*yp[1] + f*yp[2]
9.223372036854776e18
numpy clearly converts things to Float64 first, so it prevents the overflow, but this also destroys the ability to do generic computation—for example, you can’t autodiff through a computation like that. So in Julia we would definitely not want to coerce to Float64 first.
On balance I suspect that despite the slightly larger inaccuracy the Interpolations.jl approach is actually the more useful of the two.