I was playing around with Newton’s method the other day and was mind-blown when I realised that a fairly straightforward Julia implementation of Newton’s method applied to f(x) = 1/x - y
comes within roughly a factor 1.5 of inv(y)
!
"""
newton_inv(y)
Compute `inv(y)` for `1 <= y <= 2` using Newton's method.
"""
function newton_inv(y)
@fastmath x = 24/17 - (8/17) * y
for k = 1:4
@fastmath x = x * (2 - y*x)
end
return x
end
using BenchmarkTools
using Printf
function benchmark()
@printf(" inv runtime: %6.2f nanoseconds\n", @belapsed( inv($(Ref(1.0))[]))*1e9)
@printf(" Newton runtime: %6.2f nanoseconds\n", @belapsed(newton_inv($(Ref(1.0))[]))*1e9)
end
julia> BenchmarkTools.DEFAULT_PARAMETERS.samples = 1_000_000;
julia> benchmark() # Laptop plugged in
inv runtime: 1.52 nanoseconds
Newton runtime: 2.13 nanoseconds
julia> benchmark() # Laptop on battery
inv runtime: 2.53 nanoseconds
Newton runtime: 3.37 nanoseconds
Furthermore, newton_inv()
appears to produce almost correctly rounded answers:
function accuracy()
y = 1.0 .+ rand(100_000)
x = newton_inv.(y)
println(" Errors: ", unique(abs.(x .- inv.(y))))
end
julia> accuracy()
Errors: [0.0, 1.1102230246251565e-16]
Some details regarding the algorithm are described on the Wikipedia page for division algorithms.