Created
February 15, 2014 01:16
-
-
Save jwmerrill/9012954 to your computer and use it in GitHub Desktop.
Benchmark simple bisection against Roots.fzero polyalgorithm
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
using Roots | |
# Alternative "mean" definition that operates on the binary representation | |
# of a float. Using this definition, bisection will never take more than | |
# 64 steps. | |
function _middle(x::Float64, y::Float64) | |
# Use the usual float rules for combining non-finite numbers | |
if !isfinite(x) || !isfinite(y) | |
return x + y | |
end | |
# Always return 0.0 when inputs have opposite sign | |
if sign(x) != sign(y) && x != 0.0 && y != 0.0 | |
return 0.0 | |
end | |
negate = x < 0.0 || y < 0.0 | |
xint = reinterpret(Uint64, abs(x)) | |
yint = reinterpret(Uint64, abs(y)) | |
unsigned = reinterpret(Float64, (xint + yint) >> 1) | |
return negate ? -unsigned : unsigned | |
end | |
function bisect_root(fn, lower, upper) | |
x0 = lower | |
x2 = upper | |
x1 = _middle(x0, x2) | |
y0 = fn(x0) | |
y1 = fn(x1) | |
y2 = fn(x2) | |
while x0 < x1 && x1 < x2 | |
if sign(y0) == sign(y1) | |
x0, x2 = x1, x2 | |
y0, y2 = y1, y2 | |
else | |
x0, x2 = x0, x1 | |
y0, y2 = y0, y1 | |
end | |
x1 = _middle(x0, x2) | |
y1 = fn(x1) | |
end | |
return abs(y0) < abs(y2) ? x0 : x2 | |
end | |
function count_calls(method, fn, lower, upper) | |
let i = 0 | |
function g(x) | |
i += 1 | |
fn(x) | |
end | |
return method(g, lower, upper), i | |
end | |
end | |
function time_bisect() | |
accum = 0.0 | |
for i = 1:10000 | |
accum += bisect_root(sin, 3.0, 4.0) | |
end | |
return accum | |
end | |
function time_fzero() | |
accum = 0.0 | |
for i = 1:10000 | |
accum += fzero(sin, 3.0, 4.0) | |
end | |
return accum | |
end | |
# Warm up JIT | |
time_bisect() | |
time_fzero() | |
println("Number of calls for bisect_root: ", count_calls(bisect_root, sin, 3.0, 4.0)) | |
println("Number of calls for fzero: ", count_calls(fzero, sin, 3.0, 4.0)) | |
println() | |
t1 = @elapsed time_bisect() | |
t2 = @elapsed time_fzero() | |
println("Elapsed time for bisect_root: ", t1) | |
println("Elapsed time for fzero: ", t2) | |
println("Ration tfzero/tbisect_root: ", t2/t1) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
For me, this outputs:
In other words, the current algorithm uses more function evaluations, and is 8 times slower in a simple case, than bisecton.