Slow Newton method

The following code, based on @dpsanders 's hands-on Julia example seems to allocate much more than necessary, and is a lot slower as a result:

newton(f, f′, x) = x - f(x) / f′(x)

function do_newton(g, g′, x0)
    x = x0

    tolerance = 1e-15
    while abs(g(x)) > tolerance
        x = newton(g, g′, x)
    end
    
    return x
end

h(z) = z^3 - 1
h′(z) = 3z^2

function newton_fractal(width, step)
    range = -width:step:width
    L = length(range)

    N = zeros(Complex{Float64}, L, L)
    
    for I in CartesianRange(size(N))
        a = range[I[1]]
        b = range[I[2]]

        z = b + a*im

        N[I] = do_newton(h, h′, z)

        # tolerance = 1e-15
        # while abs(h(z)-0) > tolerance
        #     z = newton(h,h′,z)
        # end
        # N[I] = z
    end
    return N
end

M = newton_fractal(2.0,0.005)
@time M = newton_fractal(2.0,0.005)

Running the above code as-is results in:

1.877332 seconds (45.34 M allocations: 1.273 GiB, 1.57% gc time)

However, if I “inline” the do_newton by uncommenting the commented lines, the picture is entirely different:

0.275165 seconds (86 allocations: 9.796 MiB)

And as I typed the “inline” above, I just realized that I could put @inline in front of the do_newton function, and then both versions magically perform the same. So actually it’s solved, but I struggled too long with this not to post it, and anyway it wouldn’t have been solved without typing out the text to explain the problem. Thanks for reading, and if anyone could explain why the difference is so huge that would also be interesting.

It’s the specialization heuristic for ::Function arguments.