I have been unable to find any real code implementation of ITP method - Wikipedia. However, it claims to be the bee’s knees - just read the first paragraph of the Wikipedia article:

“In numerical analysis, the **ITP method** , short for *Interpolate Truncate and Project* , is the first root-finding algorithm that achieves the superlinear convergence of the secant method[1] while retaining the optimal[2] worst-case performance of the bisection method.[3] It is also the first method with guaranteed average performance strictly better than the bisection method under any continuous distribution.[3] In practice it performs better than traditional interpolation and hybrid based strategies (Brent’s Method, Ridders, Illinois), since it not only converges super-linearly over well behaved functions but also guarantees fast performance under ill-behaved functions where interpolations fail.[3]”

I implemented the pseudocode in the Wikipedia article, and thought I should share it. If anyone wants to put this in a package like `Roots.jl`

, `NLSolve.jl`

or `NonlinearSolve.jl`

, feel free to - I’d be thrilled.

EDIT: provided defaults, and changed argument order to make having defaults more sensible.

EDIT 2: Put the function argument first, diverging further from the Wikipedia pseudocode, but converging towards the conventional (and sensible) ordering. Also added reference to the source of hyperparameter defaults.

```
"""
ITP(f, a, b, ϵ=eps(), κ₁=0.2/(b-a), κ₂=2, n₀=1)
Find the root of the funtion `f` on the interval [a, b], given a
tolerance ϵ and hyperparameters κ₁, κ₂, n₀. The hyperparameters are
subject to the following constraints:
0 < κ₁ < Inf
1 ≤ κ₂ < 1 + (1 + √5)/2
0 < n₀ < Inf
A Julia implementation of the Interpolate, Truncate Project algorithm,
as detailed in `https://en.wikipedia.org/wiki/ITP_method`.
Hyperparameter defaults are set on the basis of the reasoning provided
in `https://docs.rs/kurbo/0.8.1/kurbo/common/fn.solve_itp.html`
"""
function ITP(f, a, b, ϵ=eps(), κ₁=0.2/(b-a), κ₂=2, n₀=1)
0 < κ₁ < Inf || error("κ₁ must be between 0 and ∞")
1 ≤ κ₂ < 1 + (1 + √5)/2 || error("κ₂ must be between 1 and 1+ (1 + √5)/2 (1 + the golden ratio)")
0 < n₀ < Inf || error("n₀ must be between 0 and ∞")
n_1div2 = ceil(Int, log2((b-a)/2ϵ))
nₘₐₓ = n_1div2 + n₀
y_a = f(a)
y_b = f(b)
sign(y_a) == sign(y_b) && error("sign(f(a)) = sign(f(b)). There is no guaranteed root in the given interval.")
j = 0
while b-a > 2ϵ
# Calculating parameters:
x_1div2 = (a+b)/2
r = ϵ*2^(nₘₐₓ - j) - (b-a)/2
δ = κ₁*(b-a)^κ₂
# Interpolation:
x_f = /(y_b*a - y_a*b, y_b-y_a)
# Truncation:
σ = sign(x_1div2 - x_f)
δ ≤ abs(x_1div2 - x_f) ? (x_t=x_f+σ*δ) : (x_t = x_1div2)
# Projection:
abs(x_t - x_1div2) ≤ r ? (x_ITP = x_t) : (x_ITP = x_1div2 - σ*r)
# Updating Interval:
y_ITP = f(x_ITP)
if y_ITP > 0
b = x_ITP
y_b = y_ITP
elseif y_ITP < 0
a = x_ITP
y_a = y_ITP
else
a = b = x_ITP
end
j += 1
end
return (a+b)/2
end
```

One could concider adding the following constraint:

` ϵ ≤ big(1e-78) && error("A tolerance this small has lead to indefinite iterations in the past. Aborting.")`