Benchmark simple bisection against Roots.fzero polyalgorithm

fzero_vs_bisect_benchmark.jl

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)

For me, this outputs:

Number of calls for bisect_root: (3.141592653589793,54)
Number of calls for fzero: (3.141592653589793,83)

Elapsed time for bisect_root: 0.170091942
Elapsed time for fzero: 1.364862163
Ration tfzero/tbisect_root: 8.02426115518159

In other words, the current algorithm uses more function evaluations, and is 8 times slower in a simple case, than bisecton.