diff --git a/.gitignore b/.gitignore
index 8d5e4ce5..1335960f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,3 +7,4 @@ test/Manifest.toml
 docs/node_modules
 *.cov
 benchmark/Manifest.toml
+.idea
diff --git a/CLAUDE.md b/CLAUDE.md
index ae756f87..398bca5b 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -55,10 +55,10 @@ ext/
 
 ### Concrete rational function types
 
-**`Barycentric{T,S}`** (barycentric.jl:13) — default method; alias `AAA = Barycentric`  
+**`Barycentric{T,S}`** (barycentric.jl:13) — alias `AAA = Barycentric`  
 Fields: `nodes`, `values`, `weights`, `w_times_f`
 
-**`Thiele{T,S}`** (thiele.jl) — continued fraction representation; alias `TCF = Thiele`  
+**`Thiele{T,S}`** (thiele.jl) — default method; continued fraction representation; alias `TCF = Thiele`  
 Fields: `nodes`, `values`, `weights`  
 Two evaluation strategies: classic and numerically stable `onediv`.
 
@@ -73,14 +73,15 @@ Fields: `polynomial::ArnoldiPolynomial`, `poles`, `residues`
 - **`ContinuumApproximation{T,S,R}`** (approximation.jl:37) — continuous domain wrapper; fields: `original`, `domain`, `fun`, `allowed`, `path`, `history`
 - **`DiscreteApproximation{T,S,R}`** (approximation.jl:75) — discrete point set wrapper; fields: `data`, `domain`, `fun`, `test_index`, `allowed`, `history`
 - **`IterationRecord{R,S,T}`** (approximation.jl:5) — convergence history entry; fields: `interpolant`, `error`, `poles`
+- **`ConvergenceStatus`** (approximation.jl:45) — why an iteration stopped; fields: `reason`, `best`, `iterations`, `error`. `reason` is one of `:converged`, `:stagnated`, `:max_degree`, `:node_failure`, `:nan_weight`, `:refinement`, `:exhausted`, `:rewound`
 
 ---
 
 ## Public API
 
 ### Approximation construction
-- `approximate(f, domain; method, max_iter, tol, allowed, refinement, stagnation)` — main entry point
-- `approximate(f, domain, poles)` — least-squares with prescribed poles
+- `approximate(f, domain, method=Thiele(); max_iter, tol, allowed, refinement, stagnation)` — main entry point; `method` is an instance (`Thiele()` default, `Barycentric()`) passed as the last positional argument
+- `approximate(f, domain, poles)` — least-squares with prescribed poles (selectable via `PartialFractions()` as the last positional argument)
 - `aaa(y, z; kwargs...)` — legacy discrete AAA (deprecated)
 
 ### Rational function queries
@@ -98,6 +99,8 @@ Fields: `polynomial::ArnoldiPolynomial`, `poles`, `residues`
 - `get_function(r)`, `domain(r)` — extract components
 - `rewind(r, index)` — revert to earlier iteration
 - `get_history(r)` — convergence history
+- `status(r)` — `ConvergenceStatus` for the run, or `nothing` if none was recorded
+- `isconverged(r)` — whether the iteration reached `tol`, as opposed to stagnating or exhausting `max_degree`
 - `test_points(r)` — test point locations
 
 ### Optimization
@@ -157,10 +160,12 @@ Fields: `polynomial::ArnoldiPolynomial`, `poles`, `residues`
 ## Design Decisions
 
 1. **Two-parameter type system**: `T` for float precision, `S` for value type — supports generic arithmetic.
-2. **Barycentric as default**: most efficient/stable; aliased `AAA` for historical compatibility.
+2. **Thiele as default**: continued-fraction method used when no selector is given; `Barycentric` (aliased `AAA`) remains available and selectable.
 3. **Continuum vs. Discrete split**: `ContinuumApproximation` and `DiscreteApproximation` reflect fundamentally different strategies.
 4. **Adaptive path discretization**: `DiscretizedPath` stores multiple refinement levels in matrix form.
 5. **`allowed` parameter**: generic function to filter pole locations; enables multiply-connected domains.
 6. **Convergence history**: optional recording enables `rewind()` and convergence plots.
+   `quitting_check` returns a `(reason, best)` tuple rather than an overloaded integer, and
+   `best_acceptable` is callable on its own so failure paths need not fake a `max_iter`.
 7. **Extension architecture**: plotting and autodiff are optional — no hard dependencies.
 8. **Precompilation workload**: uses `@compile_workload` for fast time-to-first-approximation.
diff --git a/Project.toml b/Project.toml
index 7dd8ccac..82bd9fde 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ ComplexRegions = "c64915e2-6c82-11e9-38e9-1f159a780463"
 ComplexValues = "41a84b80-6cf2-11e9-379d-9df124847946"
 GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
+IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
@@ -37,6 +38,7 @@ ComplexValues = "0.3"
 ForwardDiff = "1"
 GenericLinearAlgebra = "0.3, 0.4"
 GenericSchur = "0.5"
+IntervalSets = "0.7"
 Logging = "1"
 Makie = "0.24"
 Plots = "1"
@@ -52,10 +54,11 @@ julia = "1"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ComplexRegions = "c64915e2-6c82-11e9-38e9-1f159a780463"
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
+IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ReTest = "e0db7c4e-2690-44b9-bad6-7687da720f89"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 docs = ["CairoMakie"]
-test = ["Test", "ReTest", "LinearAlgebra", "ComplexRegions", "DoubleFloats"]
+test = ["Test", "ReTest", "LinearAlgebra", "ComplexRegions", "DoubleFloats", "IntervalSets"]
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 8433d29f..a36bba7f 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -4,15 +4,27 @@ using Logging
 # silence convergence warnings so benchmark output stays clean
 global_logger(SimpleLogger(stderr, Logging.Error))
 
+# AirspeedVelocity runs this one script against both the base revision and the PR head,
+# so it has to work with either calling convention: the old API selected the interpolant
+# type with a `method = Barycentric` keyword, the current one takes an instance
+# (`Barycentric()`) as the last positional argument. Pick the form once, at load time, so
+# no branch survives into the benchmarked call. Note the deprecation shim is deliberately
+# *not* used here: it warns on every call, which would skew the timings it appears in.
+if hasmethod(Barycentric, Tuple{})
+    approx(f, domain, M; kw...) = approximate(f, domain, M(); kw...)
+else
+    approx(f, domain, M; kw...) = approximate(f, domain; method = M, kw...)
+end
+
 const SUITE = BenchmarkGroup()
 
 # --- construction cost, by method and domain ---
 SUITE["approximate"] = BenchmarkGroup()
 for (name, method) in (("aaa", Barycentric), ("thiele", Thiele))
     g = SUITE["approximate"][name] = BenchmarkGroup()
-    g["exp_interval"] = @benchmarkable approximate(exp, $unit_interval; method = $method, allowed = true)
-    g["tanh_steep"]   = @benchmarkable approximate(x -> tanh(50x), $unit_interval; method = $method, allowed = true)
-    g["abs_circle"]   = @benchmarkable approximate(z -> abs(z - 1.0001im), $unit_circle; method = $method, allowed = true)
+    g["exp_interval"] = @benchmarkable approx(exp, $unit_interval, $method; allowed = true)
+    g["tanh_steep"]   = @benchmarkable approx(x -> tanh(50x), $unit_interval, $method; allowed = true)
+    g["abs_circle"]   = @benchmarkable approx(z -> abs(z - 1.0001im), $unit_circle, $method; allowed = true)
 end
 
 # --- construction cost on a discrete point set ---
@@ -23,19 +35,19 @@ end
 SUITE["approximate_discrete"] = BenchmarkGroup()
 for (name, method) in (("aaa", Barycentric), ("thiele", Thiele))
     g = SUITE["approximate_discrete"][name] = BenchmarkGroup()
-    g["tanh_steep"]  = @benchmarkable approximate(x -> tanh(100x), $DISCRETE_PTS; method = $method, allowed = true)
-    g["abs_shift"]   = @benchmarkable approximate(x -> abs(x + 0.5 + 0.01im), $DISCRETE_PTS; method = $method, allowed = true)
-    g["sin_recip"]   = @benchmarkable approximate(x -> sin(1 / (1.05 - x)), $DISCRETE_PTS; method = $method, allowed = true)
+    g["tanh_steep"]  = @benchmarkable approx(x -> tanh(100x), $DISCRETE_PTS, $method; allowed = true)
+    g["abs_shift"]   = @benchmarkable approx(x -> abs(x + 0.5 + 0.01im), $DISCRETE_PTS, $method; allowed = true)
+    g["sin_recip"]   = @benchmarkable approx(x -> sin(1 / (1.05 - x)), $DISCRETE_PTS, $method; allowed = true)
 end
 
 # --- evaluation cost on a fixed approximant ---
 SUITE["evaluate"] = BenchmarkGroup()
-let r = approximate(x -> tanh(50x), unit_interval, method = Barycentric, allowed = true),
+let r = approx(x -> tanh(50x), unit_interval, Barycentric; allowed = true),
     z = collect(range(-1, 1, 1000))
 
     SUITE["evaluate"]["bary_vector"] = @benchmarkable $r.($z)
 end
-let r = approximate(x -> tanh(50x), unit_interval; method = Thiele, allowed = true),
+let r = approx(x -> tanh(50x), unit_interval, Thiele; allowed = true),
     z = collect(range(-1, 1, 1000))
 
     SUITE["evaluate"]["thiele_vector"] = @benchmarkable $r.($z)
@@ -43,6 +55,6 @@ end
 
 # --- pole solve ---
 SUITE["poles"] = BenchmarkGroup()
-let r = approximate(x -> 1 / sqrt(x^2 + 0.01), unit_interval; method = Thiele)
+let r = approx(x -> 1 / sqrt(x^2 + 0.01), unit_interval, Thiele)
     SUITE["poles"]["thiele"] = @benchmarkable poles($r)
 end
diff --git a/docs/Project.toml b/docs/Project.toml
index b208be3c..2da5ca58 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -15,3 +15,4 @@ RationalFunctionApproximation = {path = ".."}
 [compat]
 Documenter = "1"
 DomainColoring = "2"
+DocumenterCitations = "~1.4"
diff --git a/docs/src/algorithms.md b/docs/src/algorithms.md
index 23b348c3..10f84e18 100644
--- a/docs/src/algorithms.md
+++ b/docs/src/algorithms.md
@@ -26,8 +26,8 @@ The AAA algorithm is the best-known and most widely used method for rational app
 The `convergenceplot` function shows the errors of the approximants found during the AAA iteration.
 
 ```@example convergence
-f = x -> cos(exp(3x))
-r = approximate(f, unit_interval)
+f(x) = cos(exp(3x))
+r = approximate(f, -1..1, AAA())
 convergenceplot(r)
 ```
 
@@ -60,8 +60,8 @@ That's why the reported error is not very large. It's worth keeping in mind that
 It is possible for the iteration to stagnate with bad poles if the original function has a singularity very close to the domain.
 
 ```@example convergence
-f = x -> tanh(500*(x - 1//4))
-r = approximate(f, unit_interval)
+f(x) = tanh(500*(x - 1//4))
+r = approximate(f, -1..1, AAA())
 convergenceplot(r)
 ```
 
@@ -69,22 +69,22 @@ This effect is thought to be mainly due to roundoff and conditioning of the prob
 
 ```@example convergence
 using DoubleFloats, ComplexRegions
-r = approximate(f, Segment{Double64}(-1, 1))
+r = approximate(f, Double64(-1)..1, AAA())
 convergenceplot(r)
 ```
 
 In the extreme case of a function with a singularity on the domain, the convergence can be substantially affected:
 
 ```@example convergence
-f = x -> abs(x - 1/8)
-r = approximate(f, unit_interval)
+f(x) = abs(x - 1/8)
+r = approximate(f, -1..1, AAA())
 convergenceplot(r)
 ```
 
 In such a case, we might get improvement by increasing the number of allowed consecutive failures via the `stagnation` keyword argument:
 
 ```@example convergence
-r = approximate(f, unit_interval, stagnation=50)
+r = approximate(f, -1..1, AAA(); stagnation=50)
 convergenceplot(r)
 ```
 
@@ -92,24 +92,24 @@ However, AAA is an $O(n^4)$ algorithm, so venturing into higher degrees can beco
 
 ## Thiele continued fractions (TCF)
 
-The TCF algorithm [SalazarCelisNumericalContinued2024](@cite) is much newer than AAA and less thoroughly battle-tested, even though it's based on a continued fraction representation of rational functions that is over a century old. Like AAA, it uses iterative greedy node selection, and the effects of that ordering look good in experiments so far but are poorly understood theoretically. In TCF's favor are its $O(n^3)$ complexity requirement and an algorithmic simplicity that requires nothing more than basic arithmetic.
+The greedy TCF algorithm [SalazarCelisNumericalContinued2024](@cite) is much newer than AAA despite being based on a continued fraction representation of rational functions that is over a century old. Like AAA, it uses iterative greedy node selection, and the effects of that ordering look very good in experiments. In TCF's favor are its $O(n^3)$ complexity requirement and an algorithmic simplicity that requires nothing more than basic arithmetic.
 
-To try greedy TCF, use `method=Thiele` or `method=TCF` as an argument to `approximate`.
+As of version 0.4 of the package, TCF is the default method used by `approximation`. You can also select it manually by passing `Thiele()` (or its alias `TCF()`) as the third positional argument.
 
 ```@example convergence
-f = x -> cos(41x - 5) * exp(-10x^2)
-r = approximate(f, unit_interval; method=TCF)
+f(x) = cos(41x - 5) * exp(-10x^2)
+r = approximate(f)    # -1..1 and TCF() by default
 convergenceplot(r)
 ```
 
-The $x$-axis of the convergence plot shows the degree of the denominator polynomial. Because the Thiele method alternates between interpolants of type $(n, n)$ and $(n+1, n)$, there are two dots in the plot for each degree. The dots corresponding to approximations of the diagonal and superdiagonal rational type are connected by lines; sometimes, they could be viewed as separate convergence curves.
+The $x$-axis of the convergence plot shows the degree of the denominator polynomial. Because the Thiele method alternates between interpolants of type $(n, n)$ and $(n+1, n)$, there are two dots in the plot for each degree, and line segments connect the dots sharing the same rational type, giving two convergence curves.
 
-Because TCF uses only addition, multiplication, and division, it is easy to use in extended precision arithmetic. Here, we use `allowed=true` to disable checking for poles, because doing so requires solving an eigenvalue problem that is far more expensive than the iteration itself.
+Because TCF uses only addition, multiplication, and division, it is easy to use in extended precision arithmetic. 
 
 ```@example convergence
-f = x -> atan(1e5*(x - 1//2))
-domain = Segment{BigFloat}(-1, 1)
-@elapsed r = approximate(f, domain; method=TCF, max_iter=400, allowed=true, stagnation=40)
+f(x) = atan(1e5*(x - 1//2))
+domain = big(-1)..1    # use BigFloats
+@elapsed r = approximate(f, domain; max_degree=150, stagnation=40)
 ```
 
 ```@example convergence
@@ -129,9 +129,9 @@ When posed on a discrete set of test points, this is a linear least-squares prob
 There is no iteration on the degree of the polynomial or rational parts of the approximant. In the continuum variant, though, the discretization of the boundary of the domain is refined iteratively until either the max-norm error is below a specified threshold or has stopped improving.
 
 ```@example convergence
-f = x -> tanh(x)
+f(x) = tanh(x)
 ζ = 1im * π * [-1/2, 1/2, -3/2, 3/2]
-r = approximate(f, Segment(-2, 2), ζ)
+r = approximate(f, -2..2, ζ)
 ```
 
 ```@example convergence
@@ -142,7 +142,7 @@ println("Max error: $(max_err(r))")
 To get greater accuracy, we can increase the degree of the polynomial part.
 
 ```@example convergence
-r = approximate(f, Segment(-2, 2), ζ; degree=20)
+r = approximate(f, -2..2, ζ; degree=20)
 max_err(r)
 ```
 
@@ -152,16 +152,16 @@ Note that the residues, which are all equal to 1 for the exact function, may not
 Pair.(residues(r)...)
 ```
 
-Suppose now we approximate $|x|$ using AAA. We can extract the poles of the result.
+Suppose now we approximate $|x|$ using TCF. We can extract the poles of the result, filtering out those that lie on the real axis.
 
 ```@example convergence
-r = approximate(abs, unit_interval, tol=1e-9)
-ζ = poles(r)
+r = approximate(abs; tol=1e-9, stagnation=30)
+ζ = filter(z -> abs(imag(z)) > 1e-8, poles(r))
 ```
 
-To what extent might these poles be suitable for a different function that has the same singularity?
+These poles might be suitable for a different function that has the same singularity:
 
 ```@example convergence
-s = approximate(x -> exp(abs(x)), unit_interval, ζ; degree=20)
+s = approximate(x -> exp(abs(x)), -1..1, ζ; degree=20)
 max_err(r)
 ```
diff --git a/docs/src/discrete.md b/docs/src/discrete.md
index d78d4149..4db9e37f 100644
--- a/docs/src/discrete.md
+++ b/docs/src/discrete.md
@@ -5,9 +5,9 @@ For many functions, discretization of the domain is straightforward. But if the
 The `approximate` function can take a vector of sample points as a domain. The given function is then evaluated only at those points, and the rational approximation is a fully discrete process that uses only the given data.
 
 ```@example mode
-using RationalFunctionApproximation, ComplexRegions
+using RationalFunctionApproximation
 x = -1:0.01:1
-f = x -> tanh(5 * (x - 0.2))
+f(x) = tanh(5 * (x - 0.2))
 r = approximate(f, x)
 ```
 You can alternatively provide just the discrete function values yourself. The domain is always given second:
@@ -20,6 +20,7 @@ r = approximate(y, x)
 As long as there are no singularities as close to the domain as the sample points are to one another, a basic discretization works well.
 
 ```@example mode
+using ComplexRegions    # to get dist() function
 I = unit_interval
 println("nearest pole is $(minimum(dist(z, I) for z in poles(r))) away")
 _, err = check(r);
@@ -29,7 +30,7 @@ println("max error on the given domain: ", maximum(abs, err))
 But if the distance to a singularity is comparable to the sample spacing, the quality of the approximation may suffer. Even worse, the method may not be aware that it has failed.
 
 ```@example mode
-f = x -> tanh(400 * (x - 0.2))
+f(x) = tanh(400 * (x - 0.2))
 r = approximate(f, x)
 println("nearest pole is $(minimum(dist(z, I) for z in poles(r))) away")
 _, err = check(r);
@@ -41,7 +42,7 @@ println("max error on finer test points: ", err)
 In the continuous mode, the adaptive sampling of the domain attempts to ensure that the approximation is accurate everywhere.
 
 ```@example mode
-r = approximate(f, I; tol=1e-12)
+r = approximate(f; tol=1e-12)
 err = maximum(abs(f(x)- r(x)) for x in range(-1, 1, 3000))
 println("max error on finer test points: ", err)
 ```
diff --git a/docs/src/domains.md b/docs/src/domains.md
index a9a10086..6087f314 100644
--- a/docs/src/domains.md
+++ b/docs/src/domains.md
@@ -10,7 +10,7 @@ The domain `unit_circle` is predefined. Here's a function approximated on the un
 using RationalFunctionApproximation, CairoMakie, DomainColoring
 const shg = current_figure
 
-f = z -> (z^3 - 1) / sin(z - 0.9 - 1im)
+f(z) = (z^3 - 1) / sin(z - 0.9 - 1im)
 r = approximate(f, unit_circle)
 ```
 
@@ -23,7 +23,6 @@ errorplot(r)
 Here is how the approximation looks in the complex plane (using a black cross to mark the pole):
 
 ```@example shapes
-using ComplexRegions
 domaincolor(r, 1.5, abs=true)
 lines!(unit_circle, color=:white, linewidth=4)
 scatter!(poles(r), markersize=16, color=:black, marker=:xcross)
@@ -36,25 +35,26 @@ Above, you can also see the zeros at roots of unity.
 This next function has infinitely many poles and an essential singularity inside the unit disk:
 
 ```@example shapes
-f = z -> tan(1 / z^4)
+f(z) = tan(1 / z^4)
 r = approximate(f, unit_circle)
 domaincolor(r, 1.5, abs=true)
 lines!(unit_circle, color=:white, linewidth=4)
 shg()
 ```
 
-We can request an approximation that is analytic in a region. In this case, it would not make sense to request one on the unit disk, since the singularities are necessary:
+We can request an approximation that is analytic in a region. The keyword `allowed=:strict` prevents the iteration from accepting a rational function with poles inside the domain. In this case, since it does not make sense to request one free of poles in the unit disk, the iteration fails:
 
 ```@example shapes
-r = approximate(f, unit_disk)
+r = approximate(f, unit_disk; allowed=:strict)
 ```
 
 In the result above, the approximation is simply a constant function, as the algorithm could do no better. However, if we request analyticity in the region exterior to the circle, everything works out:
 
 ```@example shapes
-r = approximate(f, exterior(unit_circle))
-max_err = maximum(abs, check(r, quiet=true)[2])
-println("Max error: ", max_err)
+using ComplexRegions
+r = approximate(f, exterior(unit_circle); allowed=:strict)
+maxerr = maximum(abs, check(r, quiet=true)[2])
+println("Max error: ", maxerr)
 ```
 
 ## Other shapes
@@ -63,7 +63,7 @@ We are not limited to intervals and circles! There are other shapes available in
 
 ```@example shapes
 import ComplexRegions.Shapes
-r = approximate(z -> log(0.35 + 0.4im - z), interior(Shapes.cross))
+r = approximate(z -> log(0.35 + 0.4im - z), interior(Shapes.cross); allowed=:strict)
 domaincolor(r, 1.5, abs=true)
 lines!(boundary(r.domain), color=:white, linewidth=4)
 shg()
@@ -71,7 +71,7 @@ shg()
 
 ```@example shapes
 c = Shapes.hypo(5)
-r = approximate(z -> (z+4)^(-3.5), interior(c))
+r = approximate(z -> (z+4)^(-3.5), interior(c); allowed=:strict)
 domaincolor(r, 5, abs=true)
 lines!(c, color=:white, linewidth=4)
 shg()
@@ -100,13 +100,13 @@ shg()
 It's also possible to approximate on domains with an unbounded boundary curve, but this capability is not yet automated. For example, the function
 
 ```@example shapes
-f = z -> 1 / sqrt(z - (-1 + 3im))
+f(z) = 1 / sqrt(z - (-1 + 3im))
 ```
 
 is analytic on the right half of the complex plane. In order to produce an approximation on that domain, we can transplant it to the unit disk via a Möbius transformation $\phi$:
 
 ```@example shapes
-z = cispi.(range(-1, 1, length=90))           # points on the unit circle
+z = cispi.(range(-1, 1, 90))                  # points on the unit circle
 φ = Mobius( [-1, -1im, 1], [1im, 0, -1im])    # unit circle ↦ imag axis
 extrema(real, φ.(z))
 ```
@@ -114,7 +114,7 @@ extrema(real, φ.(z))
 By composing $f$ with $\phi$, we can approximate within the disk while $f$ is evaluated only on its native domain:
 
 ```@example shapes
-r = approximate(f ∘ φ, interior(unit_circle))
+r = approximate(f ∘ φ, interior(unit_circle); allowed=:strict)
 domaincolor(r, 2, abs=true)
 lines!(unit_circle, color=:white, linewidth=4)
 scatter!(nodes(r.fun), color=:black, markersize=8)
diff --git a/docs/src/index.md b/docs/src/index.md
index 6985c408..64be24f2 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -15,18 +15,18 @@ using CairoMakie
 CairoMakie.update_theme!(size = (600, 400), fontsize=11)
 const shg = current_figure
 
-f = x -> exp(cos(4x) - sin(3x))
+f(x) = exp(cos(4x) - sin(3x))
 lines(-1..1, f)
 ```
 
-To create a rational function that approximates $f$ well on this domain, we make a call to the `approximate` function:
+To create a rational function that approximates $f$, we make a call to the `approximate` function. By default, it uses the real interval $[-1,1]$ as the domain.
 
 ```@repl interval
 using RationalFunctionApproximation
-r = approximate(f, unit_interval)
+r = approximate(f)
 ```
 
-The value of `unit_interval` is defined by the package to be the interval $[-1, 1]$. The result `r` is a type (19,19) rational approximant that can be evaluated like a function:
+The result `r` is a type (19, 18) rational approximant that can be evaluated like a function:
 
 ```@repl interval
 f(0.5) - r(0.5)
@@ -47,30 +47,16 @@ scatter!(x, 0*x, markersize = 8, color=:black)
 shg()
 ```
 
-We could choose to approximate over a wider interval:
-
-```@repl interval
-using ComplexRegions
-r = approximate(f, Segment(-2, 4))
-```
-
-Note that the degree of the rational approximant increased to capture the additional complexity.
-
 One important feature of a rational function is that it can have poles, or infinite value, at the roots of the denominator polynomial. In this case, the poles hint at where the function is most sharply peaked:
 
 ```@example interval
 poleplot(r)
 ```
 
-More typically, however, a function that is well-behaved on the real axis has a singularity structure lurking in the complex plane, and the poles of rational functions provide a unique way to cope with them. For instance, let's try approximating the hyperbolic secant function:
-
-```@example interval
-r = approximate(sech, Segment(-4, 4))
-```
-
-The sech function is smooth on the real axis but has poles on the imaginary axis at odd multiples of $i\pi/2$. The rational approximant automatically locates the poles closest to the domain:
+As another example, the sech function is smooth on the real axis but has poles on the imaginary axis at odd multiples of $i\pi/2$. An approximation over $[-4,4]$ accurately locates the poles closest to the real axis:
 
 ```@example interval
+r = approximate(sech, -4..4)
 2 * poles(r) / π
 ```
 
@@ -87,17 +73,17 @@ shg()
 A meromorphic function such as sech has only those isolated poles as singularities, and getting those right is most of the battle. By contrast, the function $\log(x + 0.05i)$ has a branch point at $x = -0.05i$ necessitating a branch cut connecting it to infinity. A rational approximant uses poles to construct a proxy branch cut:
 
 ```@example interval
-f = x -> log(x + 0.05im)
-r = approximate(f, unit_interval)
+f(x) = log(x + 0.05im)
+r = approximate(f)
 domaincolor(r, 1.2; abs=true)
 lines!(r.domain, linewidth=3, color=:white)
 shg()
 ```
 
-We close this quick introduction with approximation of $|x|$, which has a singularity on the interval. A famous result of Newman in 1964 proved that the best rational approximation of degree $n$ has root-exponential convergence.
+We close this quick introduction with approximation of $|x|$, which has a singularity on the interval. A famous result of Newman in 1964 proved that the best rational approximation of degree $n$ has root-exponential convergence. In order to get the most from the approximation, we need to tell the constructor to be stubborn about declaring the iteration stagnated.
 
 ```@example interval
-r = approximate(abs, unit_interval; tol=1e-12)
+r = approximate(abs; tol=1e-12, stagnation=50)
 convergenceplot(r)
 ```
 
@@ -105,7 +91,7 @@ We find that the nodes of the approximant are also distributed (nearly) root-exp
 
 ```@example interval
 z = filter(>(0), nodes(r))
-scatter(sort(abs.(z)), axis=(ylabel="| node |", yscale=log10,))
+scatter(sort(z), axis=(xlabel="index", ylabel="node location", xscale=sqrt, yscale=log10,))
 ```
 
 ## Feedback and contributions
diff --git a/docs/src/install.md b/docs/src/install.md
index 920f20d9..4b53d2a9 100644
--- a/docs/src/install.md
+++ b/docs/src/install.md
@@ -6,7 +6,7 @@ You can download and install this package using the general registry:
 import Pkg; Pkg.add("RationalFunctionInterpolation")
 ```
 
-If you are going to use domains other than $[-1,1]$ and the unit circle, you should also install `ComplexRegions`:
+If you are going to use domains other than real intervals and the unit circle, you should also install `ComplexRegions`:
 
 ```julia
 Pkg.add("ComplexRegions")
diff --git a/docs/src/minimax.md b/docs/src/minimax.md
index 572dcf4f..1984d561 100644
--- a/docs/src/minimax.md
+++ b/docs/src/minimax.md
@@ -7,8 +7,8 @@ For example, suppose we limit the degree of a rational interpolant of a smooth f
 ```@example minimax
 using RationalFunctionApproximation, CairoMakie
 const shg = current_figure
-f = x -> exp(cos(4x) - sin(3x))
-r = approximate(f, unit_interval, max_iter=12)
+f(x) = exp(cos(4x) - sin(3x))
+r = approximate(f, -1..1, AAA(); max_degree=12)
 ```
 
 The error varies a lot in amplitude over the interval:
@@ -27,8 +27,8 @@ errorplot(r)
 As you can see above, the error is now nearly equioscillatory over the interval. Moreover, the interpolation nodes appear to have shifted to resemble Chebyshev points of the first kind. But if we try minimax approximation on the unit circle, max-norm approximation tends to lead to equally spaced nodes:
 
 ```@example minimax
-f = z -> cos(4z) - sin(3z)
-r = approximate(f, unit_circle, max_iter=10)
+f(z) = cos(4z) - sin(3z)
+r = approximate(f, unit_circle, AAA(), max_degree=10)
 r = minimax(r, 20)
 errorplot(r, use_abs=false)
 ```
diff --git a/src/RationalFunctionApproximation.jl b/src/RationalFunctionApproximation.jl
index f22b486b..4b4d5622 100644
--- a/src/RationalFunctionApproximation.jl
+++ b/src/RationalFunctionApproximation.jl
@@ -6,6 +6,8 @@ export RFA
 using LinearAlgebra, Statistics, GenericLinearAlgebra, ComplexRegions, GenericSchur, ArnoldiVandermonde
 using PyFormattedStrings
 using PrecompileTools
+using IntervalSets: ClosedInterval, leftendpoint, rightendpoint, (..)
+export ..
 import ArnoldiVandermonde: degree
 
 export unit_interval, unit_circle, unit_disk, DiscretizedPath
@@ -18,6 +20,7 @@ export nodes, weights, degree, degrees, poles, Res, residues, roots
 include("abstract-rational.jl")
 
 export approximate, get_function, domain, check, rewind, get_history, test_points
+export ConvergenceStatus, status, isconverged
 include("approximation.jl")
 
 export Barycentric, AAA, Thiele, TCF, derivative, evaluate
@@ -38,11 +41,11 @@ include("lawson.jl")
 @setup_workload begin
     x_interval = range(-1, 1, 200)
     @compile_workload begin
-        for method in (Barycentric, Thiele)
-            approximate(sin, unit_circle; method)
+        for method in (Barycentric(), Thiele())
+            approximate(sin, unit_circle, method)
             for domain in (unit_interval, x_interval)
-                approximate(sin, domain; method)
-                approximate(cis, domain; method)
+                approximate(sin, domain, method)
+                approximate(cis, domain, method)
                 approximate(x -> 1/(x^2 + 4), domain, [2im, -2im])
             end
         end
diff --git a/src/approximation.jl b/src/approximation.jl
index 8f6aa477..f2bcc122 100644
--- a/src/approximation.jl
+++ b/src/approximation.jl
@@ -19,6 +19,45 @@ function Base.show(io::IO, ::MIME"text/plain", h::IterationRecord)
 end
 # COV_EXCL_STOP
 
+"""
+    ConvergenceStatus (type)
+
+Why an approximation iteration stopped, and which iterate it returned.
+
+# Fields
+- `reason::Symbol`: the cause of termination (see below)
+- `best::Int`: index into the history of the interpolant that was returned
+- `iterations::Int`: number of iterations completed
+- `error::Float64`: estimated error of the returned interpolant
+
+# Reasons
+- `:converged`: the error fell below `tol` with acceptable poles
+- `:stagnated`: the error plateaued for `stagnation` iterations
+- `:max_degree`: the degree budget was exhausted before reaching `tol`
+- `:node_failure`: a new node could not be added
+- `:nan_weight`: a NaN weight was encountered
+- `:refinement`: the path refinement limit was exceeded
+- `:exhausted`: all available sample values were used
+- `:rewound`: the iterate was selected by [`rewind`](@ref), not by the iteration
+
+See also [`isconverged`](@ref), [`get_history`](@ref).
+"""
+struct ConvergenceStatus
+    reason::Symbol
+    best::Int
+    iterations::Int
+    error::Float64
+end
+
+ConvergenceStatus(reason::Symbol, best::Integer, history) =
+    ConvergenceStatus(reason, best, length(history), history[best].error)
+
+# COV_EXCL_START
+function Base.show(io::IO, ::MIME"text/plain", s::ConvergenceStatus)
+    print(io, "stopped by $(s.reason) after $(s.iterations) iterations with estimated error $(round(s.error, sigdigits=4))")
+end
+# COV_EXCL_STOP
+
 abstract type AbstractApproximation{T,S} <: Function end
 
 """
@@ -33,6 +72,7 @@ Approximation of a function on a domain.
 - `allowed`: function to determine if a pole is allowed
 - `path`: a `DiscretizedPath` for the domain boundary
 - `history`: all approximations in the iteration
+- `status`: why the iteration stopped, or `nothing` if not recorded
 """
 struct ContinuumApproximation{T,S,R} <: AbstractApproximation{T,S}
     original::Function
@@ -41,6 +81,7 @@ struct ContinuumApproximation{T,S,R} <: AbstractApproximation{T,S}
     allowed::Union{Bool,Function}
     path::DiscretizedPath
     history::Union{Vector{<:IterationRecord},Nothing}
+    status::Union{ConvergenceStatus,Nothing}
 end
 
 function ContinuumApproximation(
@@ -49,15 +90,17 @@ function ContinuumApproximation(
     fun::R,
     allowed::Union{Bool,Function},
     path::DiscretizedPath,
-    history=nothing
+    history=nothing,
+    status=nothing
     ) where {T,S,R<:AbstractRationalFunction{S}}
-    return ContinuumApproximation{T,S,R}(f, domain, fun, allowed, path, history)
+    return ContinuumApproximation{T,S,R}(f, domain, fun, allowed, path, history, status)
 end
 
 (f::ContinuumApproximation)(z) = f.fun(z)
 domain(r::ContinuumApproximation) = r.domain
 get_function(r::ContinuumApproximation) = r.fun
 history(r::ContinuumApproximation) = r.history
+status(r::ContinuumApproximation) = r.status
 
 """
     DiscreteApproximation (type)
@@ -71,6 +114,7 @@ Approximation of a function on a domain.
 - `test_index`: indicator of which domain points remain as test points
 - `allowed`: function to determine if a pole is allowed
 - `history`: all approximations in the iteration
+- `status`: why the iteration stopped, or `nothing` if not recorded
 """
 struct DiscreteApproximation{T,S,R} <: AbstractApproximation{T,S}
     data::Vector{S}
@@ -79,6 +123,7 @@ struct DiscreteApproximation{T,S,R} <: AbstractApproximation{T,S}
     test_index::BitVector
     allowed::Union{Bool,Function}
     history::Union{Vector{<:IterationRecord},Nothing}
+    status::Union{ConvergenceStatus,Nothing}
     function DiscreteApproximation{T,S,R}(
         data::AbstractVector{S},
         domain::AbstractVector{T},
@@ -95,25 +140,67 @@ function DiscreteApproximation(
     fun::R,
     test_index::BitVector,
     allowed::Union{Bool,Function}=true,
-    history=nothing
+    history=nothing,
+    status=nothing
     ) where {T,S,R<:AbstractRationalFunction}
-    return DiscreteApproximation{T,float(S),typeof(fun)}(float(data), domain, fun, test_index, allowed, history)
+    return DiscreteApproximation{T,float(S),typeof(fun)}(float(data), domain, fun, test_index, allowed, history, status)
 end
 
 (f::DiscreteApproximation)(z) = f.fun(z)
 domain(r::DiscreteApproximation) = r.domain
 get_function(r::DiscreteApproximation) = r.fun
 history(r::DiscreteApproximation) = r.history
+status(r::DiscreteApproximation) = r.status
+
+"""
+    status(r::Approximation)
+
+Return the [`ConvergenceStatus`](@ref) recorded when `r` was constructed, or `nothing`
+if the construction did not record one.
+
+See also [`isconverged`](@ref).
+"""
+status
+
+"""
+    isconverged(s::ConvergenceStatus)
+    isconverged(r::Approximation)
+
+Determine whether an approximation iteration stopped because it reached the requested
+tolerance, as opposed to stagnating, exhausting its degree budget, or failing. Returns
+`false` for an approximation with no recorded status.
+
+# Examples
+```julia-repl
+julia> r = approximate(exp, unit_interval);
+
+julia> isconverged(r)
+true
+```
+
+See also [`status`](@ref), [`ConvergenceStatus`](@ref).
+"""
+isconverged(s::ConvergenceStatus) = s.reason === :converged
+isconverged(r::AbstractApproximation) = isconverged(status(r))
+isconverged(::Nothing) = false
 
 # COV_EXCL_START
 function Base.show(io::IO, ::MIME"text/plain", f::ContinuumApproximation)
     print(io, f.fun)
     print(IOContext(io, :compact=>true), " constructed on: ", f.domain)
+    _show_status(io, f.status)
 end
 
 function Base.show(io::IO, ::MIME"text/plain", f::AbstractApproximation)
     print(io, f.fun)
     print(IOContext(io, :compact=>true), " constructed from $(length(f.domain)) samples")
+    _show_status(io, f.status)
+end
+
+_show_status(io::IO, ::Nothing) = nothing
+function _show_status(io::IO, s::ConvergenceStatus)
+    print(io, "\n  ")
+    show(io, MIME"text/plain"(), s)
 end
 # COV_EXCL_STOP
 
@@ -143,23 +230,26 @@ Adaptively compute a rational interpolant on a continuous or discrete domain.
 # Arguments
 ## Continuous domain
 - `f::Function`: function to approximate
-- `domain`: curve, path, or region from ComplexRegions
+- `domain`: curve, path, or region from ComplexRegions, or an `IntervalSets.ClosedInterval`
+  (e.g. `-1..1`), which is converted to a `Segment`
 
 ## Discrete domain
 - `f::Function` or `y::AbstractVector`: function or discrete values to approximate
 - `z::AbstractVector`: domain point set
 
+## Method selection
+- `method`: instance selecting the type of rational interpolant, passed as the last positional
+  argument (`Thiele()` default, `Barycentric()`); e.g. `approximate(f, domain, Barycentric())`
+
 # Keywords
-- `method::Type`: type of rational interpolant to use (`AAA` default, `TCF`, `PartialFractions`)
-- `max_iter::Integer=150`: maximum number of iterations on node addition
+- `max_degree::Integer=100`: maximum (denominator) degree of the approximation
 - `float_type::Type`: floating point type to use for the computation¹
 - `tol::Real=1000*eps(float_type)`: relative tolerance for stopping
-- `allowed::Function`: function to determine if a pole is allowed²
+- `allowed`: no checking poles if `true`, must be outside domain if `:strict`, or use provided function
 - `refinement::Integer=3`: number of test points between adjacent nodes (continuum only)
 - `stagnation::Integer=5`: number of iterations to determine stagnation
 
 ¹Default of `float_type` is the promotion of `float(1)` and the float type of the domain.
-²Default is to disallow poles on the curve or in the interior of a continuous domain, or to accept all poles on a discrete domain. Use `allowed=true` to allow all poles.
 
 # Returns
 - `r::Approximation`: the rational interpolant
@@ -171,13 +261,13 @@ See also [`ContinuumApproximation`](@ref), [`DiscreteApproximation`](@ref), [`ch
 julia> f = x -> tanh( 40*(x - 0.15) );
 
 julia> r = approximate(f, unit_interval)
-Barycentric{Float64, Float64} rational function of type (22, 22) on the domain: Path{Float64} with 1 curve
+Thiele{Float64, Float64} rational of type (22, 22) constructed on: Segment(-1.0, 1.0)
 
 julia> ( r(0.3), f(0.3) )
-(0.9999877116508015, 0.9999877116507956)
+(0.9999877116507957, 0.9999877116507956)
 
 julia> check(r);   # accuracy over the domain
-[ Info: Max error is 1.58e-13
+[ Info: Max error is 7.78e-14
 ```
 """
 approximate(f::Function, domain)
@@ -228,58 +318,108 @@ approximate(f::Function, domain, ζ::AbstractVector)
 ##### Dispatch
 #####
 
-# Each rational type implements methods based on dispatch of Type for the first argument.
-# Here, we can call those based on a method= keyword argument instead, giving a default.
-function approximate(f::Function, domain::ComplexCurveOrPath; method::Type=Barycentric, kw...)
-     approximate(method, f, domain; kw...)
+# Deprecation shim: the pre-rename API accepted `method` as a keyword argument holding
+# a type (e.g. `method=AAA`, `method=Thiele`). The current API takes an instance as the
+# last positional argument. Extract the kwarg, warn, and return the instance + remaining
+# kwargs so callers can forward positionally.
+@noinline function _pop_deprecated_method_kw(kw)
+    m = kw[:method]
+    inst = m isa Type ? m() : m
+    Base.depwarn(
+        "Passing `method` as a keyword argument to `approximate` is deprecated. " *
+        "Pass an instance as the last positional argument instead, e.g. `approximate(f, domain, $(nameof(typeof(inst)))())`.",
+        :approximate;
+        force = true    # otherwise silent under Julia's default --depwarn=no
+    )
+    rest = Base.structdiff(NamedTuple(kw), NamedTuple{(:method,)})
+    return inst, rest
 end
 
-function approximate(y::AbstractVector, z::AbstractVector; method::Type=Barycentric, kw...)
-     approximate(method, y, z; kw...)
+# Convert an `IntervalSets.ClosedInterval` (e.g. `a..b` from IntervalSets/Makie) to a
+# `Segment`, so users can write `approximate(f, -1..1)` in place of `approximate(f, Segment(-1, 1))`.
+function approximate(f::Function, I::ClosedInterval{<:Real}, args...; kw...)
+    return approximate(f, Segment(leftendpoint(I), rightendpoint(I)), args...; kw...)
 end
 
-function approximate(
-    f::Function, domain::ComplexCurveOrPath, ζ::AbstractVector;
-    method::Type=PartialFractions,
-    kw...
-    )
-    approximate(method, f, domain, ζ; kw...)
+# Each rational type implements a method that dispatches on an instance of the type
+# (e.g. `Barycentric()`, `Thiele()`, `PartialFractions()`) passed as the last positional
+# argument. The methods here supply the default selector when none is given.
+function approximate(f::Function, domain::ComplexCurveOrPath=Segment(-1, 1); kw...)
+    if haskey(kw, :method)
+        m, rest = _pop_deprecated_method_kw(kw)
+        return approximate(f, domain, m; rest...)
+    end
+    approximate(f, domain, Thiele(); kw...)
+end
+
+function approximate(y::AbstractVector, z::AbstractVector; kw...)
+    if haskey(kw, :method)
+        m, rest = _pop_deprecated_method_kw(kw)
+        return approximate(y, z, m; rest...)
+    end
+    approximate(y, z, Thiele(); kw...)
+end
+
+function approximate(f::Function, domain::ComplexCurveOrPath, ζ::AbstractVector; kw...)
+    if haskey(kw, :method)
+        m, rest = _pop_deprecated_method_kw(kw)
+        return approximate(f, domain, ζ, m; rest...)
+    end
+    approximate(f, domain, ζ, PartialFractions(); kw...)
 end
 
 # Each rational type recognizes two signatures:
 # - f::Function, domain::ComplexCurveOrPath, [poles::AbstractVector]
 # - values::AbstractVector, test_points::AbstractVector, [poles::AbstractVector]
 
-# We fill in other convenience cases here.
+# Other convenience cases.
 
-# ::Function, ::AbstractRegion
-# Given a region as domain, we interpret poles as not being allowed in that region.
-function approximate(f::Function, R::ComplexRegions.AbstractRegion; kw...)
-    r = approximate(f, R.boundary; allowed=z->!in(z,R), kw...)
-    return ContinuumApproximation(f, R, r.fun, r.allowed, r.path, r.history)
+# ::Function, ::AbstractRegion, [selector]
+function approximate(
+    f::Function, R::ComplexRegions.AbstractRegion, method::AbstractRationalFunction=Thiele();
+    allowed=true,
+    kw...
+    )
+    if haskey(kw, :method)
+        m, rest = _pop_deprecated_method_kw(kw)
+        return approximate(f, R, m; allowed, rest...)
+    end
+    if allowed == :strict
+        # only allow poles outside the region
+        allowed = z -> !in(z,R)
+    end
+    r = approximate(f, R.boundary, method; allowed, kw...)
+    return ContinuumApproximation(f, R, r.fun, r.allowed, r.path, r.history, r.status)
 end
 
-# ::Function, ::AbstractVector
+# ::Function, ::AbstractVector, [selector]
 # Evaluate the function to call a fully discrete approximation.
 function approximate(
-    f::Function, z::AbstractVector;
+    f::Function, z::AbstractVector, method::AbstractRationalFunction=Thiele();
     allowed = true,
     kw...
     )
+    if haskey(kw, :method)
+        m, rest = _pop_deprecated_method_kw(kw)
+        return approximate(f, z, m; allowed, rest...)
+    end
     y = f.(z)
-    r = approximate(y, z; allowed, kw...)
-    return DiscreteApproximation(y, z, r.fun, r.test_index, r.allowed, r.history)
+    r = approximate(y, z, method; allowed, kw...)
+    return DiscreteApproximation(y, z, r.fun, r.test_index, r.allowed, r.history, r.status)
 end
 
-# ::Function,::AbstractVector, ::AbstractVector
+# ::Function, ::AbstractVector, ::AbstractVector, [selector]
 function approximate(
-    f::Function, z::AbstractVector, ζ::AbstractVector;
-    method = PartialFractions,
+    f::Function, z::AbstractVector, ζ::AbstractVector, method::AbstractRationalFunction=PartialFractions();
     kw...
     )
+    if haskey(kw, :method)
+        m, rest = _pop_deprecated_method_kw(kw)
+        return approximate(f, z, ζ, m; rest...)
+    end
     y = f.(z)
-    r = approximate(method, y, z, ζ; kw...)
-    return DiscreteApproximation(y, z, r.fun, r.test_index, true, r.history)
+    r = approximate(y, z, ζ, method; kw...)
+    return DiscreteApproximation(y, z, r.fun, r.test_index, true, r.history, r.status)
 end
 
 #####
@@ -317,17 +457,18 @@ Rewind a rational approximation to a state encountered during an iteration.
 # Examples
 ```jldoctest
 julia> r = approximate(x -> cos(20x), unit_interval)
-Barycentric{Float64, Float64} rational interpolant of type (24, 24) on the domain: Path{Float64} with 1 curve
+Thiele{Float64, Float64} rational of type (27, 27) constructed on: Segment(-1.0, 1.0)
 
 julia> rewind(r, 10)
-Barycentric{Float64, Float64} rational interpolant of type (10, 10) on the domain: Path{Float64} with 1 curve
+Thiele{Float64, Float64} rational of type (5, 4) constructed on: Segment(-1.0, 1.0)
 ```
 """
 function rewind(r::AbstractApproximation, idx::Integer)
     if isnothing(r.history)
         @error("No convergence history exists.")
     end
-    return typeof(r)(r.original, r.domain, r.history[idx].interpolant, r.allowed, r.path, r.history)
+    stop = ConvergenceStatus(:rewound, idx, r.history)
+    return typeof(r)(r.original, r.domain, r.history[idx].interpolant, r.allowed, r.path, r.history, stop)
 end
 
 """
@@ -420,23 +561,36 @@ function get_history(r::AbstractApproximation{T,S}; get_poles=!(r.allowed == tru
     return deg, err, zp, allowed, best
 end
 
-# Return values for quitting_check:
-#  -1: success
-#   0: continue
-#   n: iteration number to stop at
+# Index of the lowest-error iterate whose poles are all allowed. If none qualifies,
+# fall back to the final iterate. Callers that stop for their own reasons (a node that
+# could not be added, a NaN weight) use this directly to choose what to return.
+function best_acceptable(history, allowed)
+    err = [h.error for h in history]
+    if (allowed === true)
+        return argmin(i -> err[i], (i for i in eachindex(err) if !isnan(err[i])))
+    end
+    for k in sortperm(err)
+        history[k].poles = @coalesce history[k].poles poles(history[k].interpolant)
+        all(allowed, history[k].poles) && return k
+    end
+    return lastindex(err)
+end
+
+# Decide whether the iteration should stop, and on which iterate. Returns a tuple
+# `(reason, best)`; a reason of `:iterating` means carry on, and leaves `best` at 0.
+# Every other reason is one of those documented for `ConvergenceStatus`.
 function quitting_check(history, stagnation, tol, fmax, max_iter, allowed)
     n = length(history)
     err = [h.error for h in history]
 
     # Check for convergence
     # If allowed === true, do not check for allowed poles
-    status = 0
     if (err[end] <= tol*fmax)
         if (allowed === true)
-            status = -1
+            return (:converged, n)
         else
             zp = history[end].poles = poles(history[end].interpolant)
-            status = all(allowed, zp) ? -1 : 0
+            all(allowed, zp) && return (:converged, n)
         end
     end
 
@@ -448,23 +602,11 @@ function quitting_check(history, stagnation, tol, fmax, max_iter, allowed)
         stagnant = all(plateau < e for e in last(err, stagnation)) || (min_k < n - 2stagnation)
     end
 
-    # Decide on unsuccessful stopping
-    if (n >= max_iter) || stagnant
-        # Look for the best acceptable approximation:
-        if (allowed === true)
-            n = argmin(i -> err[i], (i for i in eachindex(err) if !isnan(err[i])))
-        else
-            for k in sortperm(err)
-                history[k].poles = @coalesce history[k].poles poles(history[k].interpolant)
-                if all(allowed, history[k].poles)
-                    n = k
-                    break
-                end
-            end
-        end
-        status = n
-    end
-    return status
+    # Decide on unsuccessful stopping. Stagnation is tested first, so that a run which
+    # has genuinely plateaued at the last permitted iteration is not blamed on the budget.
+    stagnant && return (:stagnated, best_acceptable(history, allowed))
+    (n >= max_iter) && return (:max_degree, best_acceptable(history, allowed))
+    return (:iterating, 0)
 end
 
 #####
@@ -479,7 +621,7 @@ Create an approximation of the derivative of `r` on the same domain.
 function derivative(r::AbstractApproximation, order=1; kwargs...)
     # TODO: This ought to be handled by dispatch on a type parameter.
     return if isa(get_function(r), AbstractRationalInterpolant)
-        approximate(derivative(get_function(r), order), domain(r); method=typeof(get_function(r)), kwargs...)
+        approximate(derivative(get_function(r), order), domain(r), get_function(r); kwargs...)
     else
         @error("Not supported. Take the derivative of the `.fun` field.")
     end
@@ -496,17 +638,17 @@ end
 
 function Base.:+(r::AbstractApproximation, g::Function)
     f(z) = r(z) + g(z)
-    return approximate(f, domain(r); method=typeof(get_function(r)))
+    return approximate(f, domain(r), get_function(r))
 end
 
 function Base.:+(r::ContinuumApproximation, s::Number)
     rs = get_function(r) + s
-    return ContinuumApproximation(rs, domain(r), rs, r.allowed, r.path, r.history)
+    return ContinuumApproximation(rs, domain(r), rs, r.allowed, r.path, r.history, r.status)
 end
 
 function Base.:+(r::DiscreteApproximation, s::Number)
     rs = get_function(r) + s
-    return DiscreteApproximation(r.data .+ s, domain(r), rs, r.test_index, r.allowed, r.history)
+    return DiscreteApproximation(r.data .+ s, domain(r), rs, r.test_index, r.allowed, r.history, r.status)
 end
 
 Base.:+(s::Union{Function,Number}, r::AbstractApproximation) = r + s
@@ -519,12 +661,12 @@ Base.:-(r::Union{Function,Number}, s::AbstractApproximation) = -s + r
 # unary -
 function Base.:-(r::ContinuumApproximation)
     rs = -get_function(r)
-    return ContinuumApproximation(rs, domain(r), rs, r.allowed, r.path, r.history)
+    return ContinuumApproximation(rs, domain(r), rs, r.allowed, r.path, r.history, r.status)
 end
 
 function Base.:-(r::DiscreteApproximation)
     rs = -get_function(r)
-    return DiscreteApproximation(-r.data, domain(r), rs, r.test_index, r.allowed, r.history)
+    return DiscreteApproximation(-r.data, domain(r), rs, r.test_index, r.allowed, r.history, r.status)
 end
 
 # * and / with 3 levels of generality
@@ -537,17 +679,17 @@ end
 
 function Base.:*(r::AbstractApproximation, g::Function)
     f(z) = r(z) * g(z)
-    return approximate(f, domain(r); method=typeof(get_function(r)))
+    return approximate(f, domain(r), get_function(r))
 end
 
 function Base.:*(r::ContinuumApproximation, s::Number)
     rs = get_function(r) * s
-    return ContinuumApproximation(rs, domain(r), rs, r.allowed, r.path, r.history)
+    return ContinuumApproximation(rs, domain(r), rs, r.allowed, r.path, r.history, r.status)
 end
 
 function Base.:*(r::DiscreteApproximation, s::Number)
     rs = get_function(r) * s
-    return DiscreteApproximation(r.data * s, domain(r), rs, r.test_index, r.allowed, r.history)
+    return DiscreteApproximation(r.data * s, domain(r), rs, r.test_index, r.allowed, r.history, r.status)
 end
 
 Base.:*(s::Union{Function,Number}, r::AbstractApproximation) = r * s
@@ -561,12 +703,12 @@ end
 
 function Base.:/(r::AbstractApproximation, g::Function)
     f(z) = r(z) / g(z)
-    return approximate(f, domain(r); method=typeof(get_function(r)))
+    return approximate(f, domain(r), get_function(r))
 end
 
 function Base.:/(r::Function, s::AbstractApproximation)
     f(z) = r(z) / s.fun(z)
-    return approximate(f, domain(s); method=typeof(s.fun))
+    return approximate(f, domain(s), s.fun)
 end
 
 Base.:/(r::AbstractApproximation, s::Number) = iszero(s) ? throw(DomainError("Division by zero")) : r * (1 / s)
@@ -575,5 +717,5 @@ Base.:/(r::Number, s::AbstractApproximation) = (z -> r) / s
 # composition
 function Base.:∘(f::Function, g::AbstractApproximation)
     # No domain checking is attempted.
-    return approximate(f ∘ g.fun, g.domain; method=typeof(g.fun))
+    return approximate(f ∘ g.fun, g.domain, g.fun)
 end
diff --git a/src/barycentric.jl b/src/barycentric.jl
index cae4ca59..ebe7ef4b 100644
--- a/src/barycentric.jl
+++ b/src/barycentric.jl
@@ -107,6 +107,9 @@ function Barycentric(points::AbstractVector, values::AbstractVector, nodeidx::Ab
     return Barycentric(points, values, idx)
 end
 
+# Empty instance, used only as a method selector in `approximate`.
+Barycentric() = Barycentric(Float64[], Float64[], Float64[])
+
 Base.copy(r::Barycentric) =
     Barycentric(copy(r.nodes), copy(r.values), copy(r.weights), copy(r.w_times_f))
 
@@ -376,19 +379,22 @@ function _initialize!(τ, fτ, C, L, f, σ, fσ, idx_test)
     return nothing
 end
 
-# TODO: This should probably enforce parameters S and T
-approximate(::Type{Barycentric{S,T}}, args...; kw...) where {S,T} = approximate(Barycentric, args...; kw...)
-
-function approximate(::Type{Barycentric},
-    f::Function, d::ComplexCurveOrPath;
+function approximate(
+    f::Function, d::ComplexCurveOrPath, ::Barycentric;
     float_type::Type = promote_type(real_type(d), typeof(float(1))),
     tol::Real = 1000*eps(float_type),
-    allowed::Union{Function,Bool} = z -> dist(z, d) > tol,
-    max_iter::Int = 150,
+    allowed = true,
+    max_degree::Int = 100,
+    max_iter = max_degree,
     refinement::Int = 3,
     stagnation::Int = 5
     )
 
+    if allowed == :strict
+        # only allow poles off the curve
+        allowed = z -> dist(z, d) > tol
+    end
+
     num_ref = 15    # initial number of test points between nodes; decreases to `refinement`
     path = DiscretizedPath(d, [0, 1]; refinement=num_ref, maxpoints=max_iter * refinement)
     σ = isclosed(d) ? points(d, [0]) : points(d, [0, 1])    # initial nodes
@@ -408,6 +414,7 @@ function approximate(::Type{Barycentric},
     # Main iteration
     idx_new_test = nothing
     n = 1       # iteration counter
+    stop = nothing
     while true
         Cmatrix = reshape(view(C, idx_test, 1:numnodes), :, numnodes)
         evaluate!(view(rτ, idx_test), r, Cmatrix)    # r at test points
@@ -415,12 +422,15 @@ function approximate(::Type{Barycentric},
         err_max, idx_max = findmax(err)
         history[n].error = err_max
 
-        status = quitting_check(history, stagnation, tol, fmax, max_iter, allowed)
-        if status > 0
-            @info("Stopping with estimated error $(round(history[status].error, sigdigits=4)) after $n iterations")
-            r = history[status].interpolant
+        reason, best = quitting_check(history, stagnation, tol, fmax, max_iter, allowed)
+        if reason !== :iterating
+            if reason !== :converged
+                @info("Stopping with estimated error $(round(history[best].error, sigdigits=4)) after $n iterations")
+                r = history[best].interpolant
+            end
+            stop = ConvergenceStatus(reason, best, history)
+            break
         end
-        (status != 0) && break
 
         ### Refinement
         idx_new = idx_test[idx_max]      # location of worst test point
@@ -429,9 +439,10 @@ function approximate(::Type{Barycentric},
             idx_new_test = add_node!(path, idx_new)
         catch
             # look for the best acceptable case
-            status = quitting_check(history, stagnation, tol, fmax, 1, allowed)
-            r = history[status].interpolant
-            @info("Unable to add new node; stopping with estimated error $(round(history[status].error, sigdigits=4))")
+            best = best_acceptable(history, allowed)
+            r = history[best].interpolant
+            stop = ConvergenceStatus(:node_failure, best, history)
+            @info("Unable to add new node; stopping with estimated error $(round(history[best].error, sigdigits=4))")
             break
         end
 
@@ -452,15 +463,16 @@ function approximate(::Type{Barycentric},
         n += 1
         numnodes += 1
     end
-    return ContinuumApproximation(f, d, r, allowed, path, history)
+    return ContinuumApproximation(f, d, r, allowed, path, history, stop)
 end
 
-function approximate(::Type{Barycentric},
-    y::AbstractVector{T}, z::AbstractVector{S};
+function approximate(
+    y::AbstractVector{T}, z::AbstractVector{S}, ::Barycentric;
     float_type::Type = promote_type(real_type(eltype(z)), typeof(float(1))),
     tol::AbstractFloat = 1000*eps(float_type),
     allowed::Union{Function,Bool} = true,
-    max_iter::Int = 100,
+    max_degree::Int = min(length(y), 100),
+    max_iter = max_degree,
     stagnation::Int = 5,
     ) where {T<:Number,S<:Number}
 
@@ -482,6 +494,7 @@ function approximate(::Type{Barycentric},
     r = Barycentric([z[i₀]], [y[i₀]], view(L, idx_test, 1:1))
     history = [IterationRecord(r, NaN, missing)]
     n = 1    # iteration counter
+    stop = nothing
     while count(idx_test) > 0
         evaluate!(view(values, idx_test), r, view(C, idx_test, 1:n))    # r at test points
         idx_max, err_max = 0, -Inf
@@ -496,12 +509,15 @@ function approximate(::Type{Barycentric},
         end
         history[n].error = err_max
 
-        status = quitting_check(history, stagnation, tol, fmax, max_iter, allowed)
-        if status > 0
-            @info("Stopping with estimated error $(round(history[status].error, sigdigits=4)) after $n iterations")
-            r = history[status].interpolant
+        reason, best = quitting_check(history, stagnation, tol, fmax, max_iter, allowed)
+        if reason !== :iterating
+            if reason !== :converged
+                @info("Stopping with estimated error $(round(history[best].error, sigdigits=4)) after $n iterations")
+                r = history[best].interpolant
+            end
+            stop = ConvergenceStatus(reason, best, history)
+            break
         end
-        (status != 0) && break
 
         # Add new node:
         idx_test[idx_max] = false
@@ -509,7 +525,8 @@ function approximate(::Type{Barycentric},
         push!(history, IterationRecord(r, NaN, missing))
         n += 1
     end
-    return DiscreteApproximation(y, z, r, idx_test, allowed, history)
+    stop = @something stop ConvergenceStatus(:exhausted, lastindex(history), history)
+    return DiscreteApproximation(y, z, r, idx_test, allowed, history, stop)
 end
 
 # Operations with scalars that can be done quickly.
diff --git a/src/lawson.jl b/src/lawson.jl
index 462253ce..41a949e3 100644
--- a/src/lawson.jl
+++ b/src/lawson.jl
@@ -103,7 +103,7 @@ function brasil(r::ContinuumApproximation{T}; tol=1000*eps(T), σmax=0.1, τ=0.1
         ℓ = c .* diff(t)
         t = [0; cumsum(ℓ) / sum(ℓ)]
         z = point(p, t[idx])
-        fun, _ = approximate(r.original.(z), z; method=Thiele, allowed=r.allowed)
+        fun, _ = approximate(r.original.(z), z, Thiele(); allowed=r.allowed)
         @show maximum(δ) / minimum(δ) - 1, extrema(δ)
         iter += 1
     end
diff --git a/src/logo.jl b/src/logo.jl
index ef076fe7..3aaf5521 100644
--- a/src/logo.jl
+++ b/src/logo.jl
@@ -1,7 +1,7 @@
 @usingany CairoMakie, Colors
 ff(x) = sqrt(x^2-1)
 # r = approximate(z->log(z+0.05+0.05im), unit_interval)
-r = approximate(ff, Segment{Double64}(101//100+1/20im,3//2+1/20im); method=Thiele, allowed=true, tol=1e-30)
+r = approximate(ff, Segment{Double64}(101//100+1/20im,3//2+1/20im), Thiele(); allowed=true, tol=1e-30)
 jgreen = RGB(0.22, 0.596, 0.149)
 jred = RGB(0.796,0.235,0.2)
 jpurple = RGB(0.584, 0.345, 0.698)
diff --git a/src/parfrac.jl b/src/parfrac.jl
index 9668dcc0..9139520d 100644
--- a/src/parfrac.jl
+++ b/src/parfrac.jl
@@ -20,7 +20,7 @@ struct PartialFractions{S} <: AbstractRationalFunction{S}
 end
 
 function PartialFractions(
-    p::ArnoldiPolynomial = ArnoldiPolynomial(),
+    p::ArnoldiPolynomial = ArnoldiPolynomial([0], ArnoldiBasis([1], 0)),
     poles::AbstractVector = ComplexF64[],
     residues::AbstractVector = ComplexF64[]
     )
@@ -111,8 +111,8 @@ function refine_by_singularity(d::ComplexCurveOrPath, ζ::AbstractVector;
     return path
 end
 
-function approximate(::Type{PartialFractions},
-    f::Function, d::ComplexCurveOrPath, ζ::AbstractVector;
+function approximate(
+    f::Function, d::ComplexCurveOrPath, ζ::AbstractVector, ::PartialFractions;
     degree = max(1, div(length(ζ), 2)),
     init =  max(400, length(d) * 100),
     refinement = 3,
@@ -125,8 +125,8 @@ function approximate(::Type{PartialFractions},
     return ContinuumApproximation(f, d, r, true, path, nothing)
 end
 
-function approximate(::Type{PartialFractions},
-    y::AbstractVector, z::AbstractVector, ζ::AbstractVector;
+function approximate(
+    y::AbstractVector, z::AbstractVector, ζ::AbstractVector, ::PartialFractions;
     degree = max(1, div(length(ζ), 2)),
     )
     r = PartialFractions(z, y, ζ, degree)
diff --git a/src/thiele.jl b/src/thiele.jl
index ddd73e5f..92d13011 100644
--- a/src/thiele.jl
+++ b/src/thiele.jl
@@ -111,6 +111,9 @@ See also [`set_eval_method`](@ref).
 """
 set_weight_method(m::ThieleMethod) = (@eval default_weight_method() = $m; m)
 
+# Empty instance, used only as a method selector in `approximate`.
+Thiele() = Thiele(Float64[], Float64[], Float64[])
+
 # Evaluation at a point
 function evaluate(r::Thiele, z::Number, method::ThieleMethod=default_eval_method())
     return if isinf(z)
@@ -435,19 +438,22 @@ function _sweep!(fτ, f::Function, τ, idx)
     return fmax
 end
 
-# TODO: This should probably enforce parameters S and T
-approximate(::Type{Thiele{S,T}}, args...; kw...) where {S,T} = approximate(Thiele, args...; kw...)
-
-function approximate(::Type{Thiele},
-    f::Function, d::Union{ComplexPath,ComplexCurve};
+function approximate(
+    f::Function, d::Union{ComplexPath,ComplexCurve}, ::Thiele;
     float_type::Type = promote_type(real_type(d), typeof(float(1))),
     tol::Real = 1000*eps(float_type),
-    allowed::Union{Function,Bool} = z -> dist(z, d) > tol,
-    max_iter::Int = 240,
+    allowed = true,
+    max_degree = 100,
+    max_iter::Int = 2max_degree,
     refinement::Int = 3,
     stagnation::Int = 5
     )
 
+    if allowed == :strict
+        # only allow poles off the curve
+        allowed = z -> dist(z, d) > tol
+    end
+
     num_ref = 15    # initial number of test points between nodes; decreases to `refinement`
     path = DiscretizedPath(d, [0, 1]; refinement=num_ref, maxpoints=max_iter * refinement)
     σ = [point(d, 0)]
@@ -491,16 +497,20 @@ function approximate(::Type{Thiele},
 
     # Main iteration
     n = 1       # iteration counter
+    stop = nothing
     while true
         err_max, idx_max = _sweep!(rbuf, zbuf, abuf, bbuf, r, τ, fτ, active)
         history[n].error = err_max
 
-        status = quitting_check(history, stagnation, tol, fmax, max_iter, allowed)
-        if status > 0
-            @info("Stopping with estimated error $(round(history[status].error, sigdigits=4)) after $n iterations")
-            r = history[status].interpolant
+        reason, best = quitting_check(history, stagnation, tol, fmax, max_iter, allowed)
+        if reason !== :iterating
+            if reason !== :converged
+                @info("Stopping with estimated error $(round(history[best].error, sigdigits=4)) after $n iterations")
+                r = history[best].interpolant
+            end
+            stop = ConvergenceStatus(reason, best, history)
+            break
         end
-        (status != 0) && break
 
         # Add node to approximant
         idx_new = active[idx_max]      # location of worst test point
@@ -509,9 +519,10 @@ function approximate(::Type{Thiele},
             push!(history, IterationRecord(r, NaN, missing))
         catch(e)
             # look for the best acceptable case
-            status = quitting_check(history, stagnation, tol, fmax, 1, allowed)
-            r = history[status].interpolant
-            @info("NaN weight encountered; stopping with estimated error $(round(history[status].error, sigdigits=4))")
+            best = best_acceptable(history, allowed)
+            r = history[best].interpolant
+            stop = ConvergenceStatus(:nan_weight, best, history)
+            @info("NaN weight encountered; stopping with estimated error $(round(history[best].error, sigdigits=4))")
             @debug("Error $e")
             break
         end
@@ -522,9 +533,10 @@ function approximate(::Type{Thiele},
             idx_new_test = add_node!(path, CartesianIndices(τ)[idx_new])
         catch
             # look for the best acceptable case
-            status = quitting_check(history, stagnation, tol, fmax, 1, allowed)
-            r = history[status].interpolant
-            @info("Maximum path refinement exceeded; stopping with estimated error $(round(history[status].error, sigdigits=4))")
+            best = best_acceptable(history, allowed)
+            r = history[best].interpolant
+            stop = ConvergenceStatus(:refinement, best, history)
+            @info("Maximum path refinement exceeded; stopping with estimated error $(round(history[best].error, sigdigits=4))")
             break
         end
 
@@ -548,15 +560,15 @@ function approximate(::Type{Thiele},
             end
         end
     end
-    return ContinuumApproximation(f, d, r, allowed, path, history)
+    return ContinuumApproximation(f, d, r, allowed, path, history, stop)
 end
 
-function approximate(::Type{Thiele},
-    y::AbstractVector{T}, z::AbstractVector{S};
+function approximate(
+    y::AbstractVector{T}, z::AbstractVector{S}, ::Thiele;
     float_type::Type = promote_type(real_type(eltype(z)), typeof(float(1))),
     tol::AbstractFloat = 1000*eps(float_type),
     allowed::Union{Function,Bool} = true,
-    max_iter::Int = length(y),
+    max_iter::Int = min(length(y), 200),
     stagnation::Int = 5,
     ) where {T<:Number,S<:Number}
 
@@ -573,6 +585,7 @@ function approximate(::Type{Thiele},
 
     history = [IterationRecord(r, NaN, missing)]
     n = 1    # iteration counter
+    stop = nothing
     while length(z) > 0
         evaluate!(r_test, r, z_test)
         @inbounds for i in eachindex(r_test)   # array evaluation skips the underflow check
@@ -582,17 +595,21 @@ function approximate(::Type{Thiele},
         err_max, idx_max = findmax(abs(e) for e in r_test)
         history[n].error = err_max
 
-        status = quitting_check(history, stagnation, tol, fmax, max_iter, allowed)
-        if status > 0
-            if isinf(err_max)
-                @info("Used all sample values without convergence")
-                status = max_iter
-            else
-                @info("Stopping with estimated error $(round(history[status].error, sigdigits=4)) after $n iterations")
+        reason, best = quitting_check(history, stagnation, tol, fmax, max_iter, allowed)
+        if reason !== :iterating
+            if reason !== :converged
+                # An infinite error estimate means there is nothing left to test against.
+                if isinf(err_max)
+                    reason, best = :exhausted, lastindex(history)
+                    @info("Used all sample values without convergence")
+                else
+                    @info("Stopping with estimated error $(round(history[best].error, sigdigits=4)) after $n iterations")
+                end
+                r = history[best].interpolant
             end
-            r = history[status].interpolant
+            stop = ConvergenceStatus(reason, best, history)
+            break
         end
-        (status != 0) && break
 
         # Add new node:
         try
@@ -603,15 +620,16 @@ function approximate(::Type{Thiele},
             deleteat!(r_test, idx_max)
         catch(e)
             # look for the best acceptable case
-            status = quitting_check(history, stagnation, tol, fmax, 1, allowed)
-            r = history[status].interpolant
-            @info("Adding node failed; stopping with estimated error $(round(history[status].error, sigdigits=4))")
+            best = best_acceptable(history, allowed)
+            r = history[best].interpolant
+            stop = ConvergenceStatus(:node_failure, best, history)
+            @info("Adding node failed; stopping with estimated error $(round(history[best].error, sigdigits=4))")
             @debug("Error $e")
             break
         end
         n += 1
     end
-    return DiscreteApproximation(y, z, r, idx_test, allowed, history)
+    return DiscreteApproximation(y, z, r, idx_test, allowed, history, stop)
 end
 
 # Operations with scalars that can be done quickly.
diff --git a/test/RFATests.jl b/test/RFATests.jl
index bb35f18e..28c07f6e 100644
--- a/test/RFATests.jl
+++ b/test/RFATests.jl
@@ -1,6 +1,6 @@
 module RFATests
 
-using RationalFunctionApproximation, ReTest, ComplexRegions, DoubleFloats, Logging
+using RationalFunctionApproximation, ReTest, ComplexRegions, DoubleFloats, IntervalSets, Logging
 const RFA = RationalFunctionApproximation
 
 pass(f, r, z; kw...) = isapprox(f.(z), r.(z), norm=u->maximum(abs, u); kw...)
@@ -16,5 +16,6 @@ include("circle.jl")
 include("custom.jl")
 include("operations.jl")
 include("parfrac.jl")
+include("deprecated.jl")
 
 end
diff --git a/test/circle.jl b/test/circle.jl
index 5bf9905c..153b4e9d 100644
--- a/test/circle.jl
+++ b/test/circle.jl
@@ -4,7 +4,7 @@
     UC = unit_circle
 
     @testset "Unit disk for $method" for method in (Barycentric, Thiele)
-        approx(f; kw...) = approximate(f, UD; method, kw...)
+        approx(f; kw...) = approximate(f, UD, method(); kw...)
         f = z -> sin(10z) * exp(-z^2); @test pass(f, approx(f), pts, rtol=2e-11)
         f = z -> sin(1/(1.1 - z)); @test pass(f, approx(f), pts, rtol=2e-13)
         f = sec; @test pass(f, approx(f, max_iter=15), pts, rtol=1e-6)
@@ -13,23 +13,23 @@
     end
 
     @testset "Unit circle for $method" for method in (Barycentric, Thiele)
-        f = z -> abs(z-1im);  @test pass(f, approximate(f, UC), pts, rtol=2e-10)
-        f = z -> tan(π*z);  @test pass(f, approximate(f, UC), pts, rtol=2e-13)
-        f = z -> tanh(100z); @test pass(f, approximate(f, UC), pts, rtol=2e-13)
+        f = z -> abs(z-1im);  @test pass(f, approximate(f, UC, Barycentric()), pts, rtol=2e-10)
+        f = z -> tan(π*z);  @test pass(f, approximate(f, UC, Barycentric()), pts, rtol=2e-13)
+        f = z -> tanh(100z); @test pass(f, approximate(f, UC, Barycentric()), pts, rtol=2e-13)
     end
 
     @testset "Array evaluation for Thiele" begin
         f = z -> sin(10z) * exp(-z^2)
-        r = approximate(f, UC, method=Thiele)
+        r = approximate(f, UC, Thiele())
         @test isapprox(f.(pts), r(pts), norm=u->maximum(abs, u), rtol=2e-11)
         f = z -> real(tan(π*z))
-        r = approximate(f, UC, method=Thiele)
+        r = approximate(f, UC, Thiele())
         @test isapprox(f.(pts), r(pts), norm=u->maximum(abs, u), rtol=2e-11)
     end
 
     @testset "Float type conversion for $method" for method in (Barycentric, Thiele)
         f = z -> abs(z - 1im)
-        r = approximate(f, UC; method)
+        r = approximate(f, UC, method())
         r32 = convert(Float32, r.fun)
         @test r32 isa method{Float32,ComplexF32}
     end
@@ -37,47 +37,47 @@
     @testset "Translate and scale for $method" for method in (Barycentric, Thiele)
         f = z -> sin(10z) * exp(-z^2)
         for (a, c) in ( (2.5, 0), (1, -1im), (0.4, -2))
-            F = approximate(f, a*UC + c)
+            F = approximate(f, a*UC + c, Barycentric())
             @test pass(f, F, a*pts .+ c)
         end
-        f = z -> 1e100sin(z); @test pass(f, approximate(f, UD), pts, rtol=2e-13)
-        @test pass(f, approximate(f, UD, max_iter=12), pts, rtol=1e-6)
+        f = z -> 1e100sin(z); @test pass(f, approximate(f, UD, Barycentric()), pts, rtol=2e-13)
+        @test pass(f, approximate(f, UD, Barycentric(); max_iter=12), pts, rtol=1e-6)
     end
 
     @testset "Poles, zeros, residues in $T for $method" for T in (Float64, Double64), method in (Barycentric, Thiele)
         UC = Circle{T}(0, 1)
         UD = interior(UC)
-        f = z -> tan(T(π)*z);  F = approximate(f, UC)
+        f = z -> tan(T(π)*z);  F = approximate(f, UC, Barycentric())
         pol = poles(F); @test sort(abs.(pol))[1:5] ≈ 0.5*[1;1;3;3;5] atol=1e-3
 
-        f = z -> exp(exp(z)) / (z - 1im // 5); pol = poles(approximate(f, UC));
+        f = z -> exp(exp(z)) / (z - 1im // 5); pol = poles(approximate(f, UC, Barycentric()));
         @test minimum(@. abs(pol - 1im // 5)) < 1000eps(T)
 
-        f = z -> (z+1) * (z+2) / ((z+3) * (z+4));  F = approximate(f, UC)
+        f = z -> (z+1) * (z+2) / ((z+3) * (z+4));  F = approximate(f, UC, Barycentric())
         pol = poles(F);  zer = roots(F);
         @test isapprox(sum(pol+zer), -10, atol=1000eps(T))
 
-        f = z -> 2/(3+z) + 5im / (z-2im);  F = approximate(f, UD)
+        f = z -> 2/(3+z) + 5im / (z-2im);  F = approximate(f, UD, Barycentric())
         @test isapprox( prod(residues(F)[2]), 10im, atol=sqrt(eps(T)) )
 
-        f = z -> (z-(3+3im))/(z+2);  F = approximate(f, UD)
+        f = z -> (z-(3+3im))/(z+2);  F = approximate(f, UD, Barycentric())
         pol, zer = poles(F), roots(F);  @test isapprox(pol[1]*zer[1], -6-6im, atol=1000eps(T))
     end
 
     @testset "Tolerance" begin
-        f = z -> exp(3*z); @test !pass(f, approximate(f, UD, tol=1e-4), pts, atol=1e-8)
-        f = z -> exp(3*z); @test pass(f, approximate(f, UD, tol=1e-10), pts, atol=1e-8)
+        f = z -> exp(3*z); @test !pass(f, approximate(f, UD, Barycentric(); tol=1e-4), pts, atol=1e-8)
+        f = z -> exp(3*z); @test pass(f, approximate(f, UD, Barycentric(); tol=1e-10), pts, atol=1e-8)
     end
 
     @testset "Low degree" begin
-        f = x -> 0; @test pass(f, approximate(f, max_iter=1, UD), pts, atol=2e-13)
-        f = x -> x; @test pass(f, approximate(f, max_iter=2, UD), pts, atol=2e-13)
-        f = x -> x+x^2; @test pass(f, approximate(f, max_iter=3, UD), pts, atol=2e-13)
-        f = x -> x+x^3; @test pass(f, approximate(f, max_iter=4, UD), pts, atol=2e-13)
-        f = x -> x+x^3; @test !pass(f, approximate(f, max_iter=3, UD), pts, atol=2e-13)
-        f = x -> 1/(3im + x + x^2); @test pass(f, approximate(f, max_iter=3, UC), pts, rtol=2e-13)
-        f = x -> 1/(3im + x + x^2); @test !pass(f, approximate(f, max_iter=2, UC), pts, rtol=2e-13)
-        f = x -> 1/(1.01 + x^3); @test pass(f, approximate(f, max_iter=4, UD), pts, rtol=2e-13)
-        f = x -> 1/(1.01 + x^3); @test !pass(f, approximate(f, max_iter=3, UD), pts, rtol=2e-13)
+        f = x -> 0; @test pass(f, approximate(f, UD, Barycentric(); max_iter=1), pts, atol=2e-13)
+        f = x -> x; @test pass(f, approximate(f, UD, Barycentric(); max_iter=2), pts, atol=2e-13)
+        f = x -> x+x^2; @test pass(f, approximate(f, UD, Barycentric(); max_iter=3), pts, atol=2e-13)
+        f = x -> x+x^3; @test pass(f, approximate(f, UD, Barycentric(); max_iter=4), pts, atol=2e-13)
+        f = x -> x+x^3; @test !pass(f, approximate(f, UD, Barycentric(); max_iter=3), pts, atol=2e-13)
+        f = x -> 1/(3im + x + x^2); @test pass(f, approximate(f, UC, Barycentric(); max_iter=3), pts, rtol=2e-13)
+        f = x -> 1/(3im + x + x^2); @test !pass(f, approximate(f, UC, Barycentric(); max_iter=2), pts, rtol=2e-13)
+        f = x -> 1/(1.01 + x^3); @test pass(f, approximate(f, UD, Barycentric(); max_iter=4), pts, rtol=2e-13)
+        f = x -> 1/(1.01 + x^3); @test !pass(f, approximate(f, UD, Barycentric(); max_iter=3), pts, rtol=2e-13)
     end
 end
diff --git a/test/deprecated.jl b/test/deprecated.jl
new file mode 100644
index 00000000..60daf4dd
--- /dev/null
+++ b/test/deprecated.jl
@@ -0,0 +1,35 @@
+# The pre-rename API selected the interpolant type with a `method` keyword argument.
+# That spelling still works, but must announce itself: `Base.depwarn` is silent under
+# Julia's default `--depwarn=no`, so the shim passes `force=true`.
+@testset "Deprecated `method` keyword" begin
+    f = exp
+    z = collect(range(-1, 1, 200))
+    warned = (:warn, r"deprecated")
+
+    @testset "warns on every entry point" begin
+        @test_logs warned match_mode=:any approximate(f, Segment(-1, 1); method=Barycentric)
+        @test_logs warned match_mode=:any approximate(f, -1..1; method=Barycentric)
+        @test_logs warned match_mode=:any approximate(f, interior(Circle(0, 1)); method=Barycentric)
+        @test_logs warned match_mode=:any approximate(f, z; method=Barycentric)
+        @test_logs warned match_mode=:any approximate(f.(z), z; method=Barycentric)
+    end
+
+    @testset "selects the requested type" begin
+        @test approximate(f, Segment(-1, 1); method=Barycentric).fun isa Barycentric
+        @test approximate(f, Segment(-1, 1); method=Thiele).fun isa Thiele
+        # an instance is accepted as well as a type
+        @test approximate(f, Segment(-1, 1); method=Barycentric()).fun isa Barycentric
+    end
+
+    @testset "other keywords survive the shim" begin
+        r = approximate(x -> exp(3x), Segment(-1, 1); method=Barycentric, tol=1e-5)
+        @test r.fun isa Barycentric
+        @test !pass(x -> exp(3x), r, range(-1, 1, 500), atol=1e-7)
+        @test pass(x -> exp(3x), r, range(-1, 1, 500), atol=5e-5)
+    end
+
+    @testset "no warning without the keyword" begin
+        @test_logs approximate(f, Segment(-1, 1))
+        @test_logs approximate(f, Segment(-1, 1), Barycentric())
+    end
+end
diff --git a/test/discrete.jl b/test/discrete.jl
index 0836658c..fd5b64de 100644
--- a/test/discrete.jl
+++ b/test/discrete.jl
@@ -4,7 +4,7 @@
         tol = 2000*eps(T)
         z = T(10) .^ range(T(-15), T(0), 500);
         pts = [-reverse(z); 0; z]
-        approx(f; kw...) = approximate(f, pts; method, kw...)
+        approx(f; kw...) = approximate(f, pts, method(); kw...)
         f = x -> abs(x + 1//2 + 1im//100); @test pass(f, approx(f), pts; rtol=tol)
         f = x -> sin(1 / (21//20 - x)); @test pass(f, approx(f), pts; rtol=tol)
         f = x -> 1im*x + exp(-1 / x^2); @test pass(f, approx(f), pts; rtol=tol)
@@ -26,7 +26,7 @@
         tol = 2000*eps(T)
         z = T(10) .^ range(T(-15), T(0), 500);
         pts = [-reverse(z); 0; z]
-        approx(f; kw...) = approximate(f, pts; method, kw...)
+        approx(f; kw...) = approximate(f, pts, method(); kw...)
         f = x -> sin(1 / (21//20 - x)); @test pass(f, approx(f), pts; rtol=tol)
         f = x -> 1im*x + exp(-1 / x^2); @test pass(f, approx(f; stagnation=30), pts; rtol=tol)
         f = x -> x + sin(80x) * exp(-10x^2); @test pass(f, approx(f; stagnation=30), pts; rtol=tol)
@@ -37,7 +37,7 @@
 
     @testset "Discrete circle for $method" for method in (Barycentric, Thiele)
         pts = cispi.(2 * (0:999) / 1000)
-        approx(f; kw...) = approximate(f, pts; method, kw...)
+        approx(f; kw...) = approximate(f, pts, method(); kw...)
         f = z -> sin(10z) * exp(-z^2); @test pass(f, approx(f), pts, rtol=2e-11)
         f = z -> sin(1/(1.1 - z)); @test pass(f, approx(f), pts, rtol=2e-13)
         f = sec; @test pass(f, approx(f), pts, rtol=1e-6)
@@ -50,7 +50,7 @@
 
     @testset "Discrete interval, low accuracy" begin
         pts = range(-1, 1, 1001)
-        approx(f; kw...) = approximate(f, pts; method=Barycentric, kw...)
+        approx(f; kw...) = approximate(f, pts, Barycentric(); kw...)
         f = x -> exp(3x);
         r = approx(f, tol=1e-5)
         @test !pass(f, r, pts, atol=1e-9)
@@ -59,7 +59,7 @@
     end
 
     @testset "Poles, zeros, residues in $T" for T in (Float64,Double64)
-        approx(f; kw...) = approximate(f, range(T(-1), T(1), 1001); method=Barycentric, kw...)
+        approx(f; kw...) = approximate(f, range(T(-1), T(1), 1001), Barycentric(); kw...)
         f = z -> (z+1) * (z+2) / ((z+3) * (z+4))
         r = approx(f)
         pol = poles(r)
@@ -81,7 +81,7 @@
 
     @testset "Vertical scaling in $T" for T in (Float64, Double64)
         pts = range(T(-1), T(1), 1001)
-        approx(f; kw...) = approximate(f, pts; method=Barycentric, kw...)
+        approx(f; kw...) = approximate(f, pts, Barycentric(); kw...)
         f = x -> T(10)^50 * sin(x); @test pass(f, approx(f), pts, rtol=2000*eps(T))
         f = x -> T(10)^(-50) * cos(x); @test pass(f, approx(f), pts, rtol=2000*eps(T))
     end
@@ -90,7 +90,7 @@
     @testset "Polynomials and reciprocals" begin
         pts = range(-1, 1, 1001)
         tol = 2000*eps(Float64)
-        approx(f; kw...) = approximate(f, pts; method=Barycentric, kw...)
+        approx(f; kw...) = approximate(f, pts, Barycentric(); kw...)
         f = x -> 0; @test pass(f, approx(f), pts, atol=tol)
         f = x -> x; @test pass(f, approx(f), pts, atol=tol)
         f = x -> 1im*x; @test pass(f, approx(f), pts, atol=tol)
@@ -105,7 +105,7 @@
     @testset "Limited degree" begin
         tol = 2000*eps(Float64)
         pts = range(-1, 1, 1001)
-        approx(f; kw...) = approximate(f, pts; method=Barycentric, kw...)
+        approx(f; kw...) = approximate(f, pts, Barycentric(); kw...)
         f = x -> 0; @test pass(f, approx(f, max_iter=1), pts, atol=tol)
         f = x -> x; @test pass(f, approx(f, max_iter=2), pts, atol=tol)
         f = x -> 1im*x; @test pass(f, approx(f, max_iter=4), pts, atol=tol)
@@ -123,7 +123,7 @@
 
     @testset "Interval [$a, $b]" for (a, b) in ((-2, 3), (0, 4), (-2e-4, 0), (-3e3, 5e6))
         pts = range(a, b, 1000)
-        approx(f; kw...) = approximate(f, pts; method=Barycentric, kw...)
+        approx(f; kw...) = approximate(f, pts, Barycentric(); kw...)
         toler(f) = 1000 * eps() * max(b - a, abs(f(a)), abs(f(b)))
         f = x -> 1 / sin(a + (b-a)*1.05im - x); @test pass(f, approx(f), pts, atol=toler(f))
         f = x -> exp(-10/(a + 1.1*(b-a) - x)); @test pass(f, approx(f), pts, atol=toler(f))
diff --git a/test/imag_interval.jl b/test/imag_interval.jl
index 65e567e6..e770f1fd 100644
--- a/test/imag_interval.jl
+++ b/test/imag_interval.jl
@@ -13,7 +13,7 @@
         T = Float64
         tol = 8000eps(T)
         pts = test_points[T]
-        approx(f; kw...) = approximate(f, domain[T]; method, kw...)
+        approx(f; kw...) = approximate(f, domain[T], method(); kw...)
         @testset "Function $iter" for (iter, f) in enumerate((
             x -> abs(x - 1//2 + 1im//100),
             x -> sinh(1 / (21//20 - x)),
@@ -34,7 +34,7 @@
         T = Double64
         pts = test_points[T]
         tol = 3000*eps(T)
-        approx(f; kw...) = approximate(f, domain[T]; method, kw...)
+        approx(f; kw...) = approximate(f, domain[T], method(); kw...)
         @testset "Function $iter" for (iter, f) in enumerate((
             x -> abs(x - 1//2 + 1im//100),
             x -> sinh(1 / (21//20 - x)),
@@ -55,7 +55,7 @@
         T = Float64
         method = Barycentric
         pts = test_points[T]
-        approx(f; kw...) = approximate(f, domain[T]; method, kw...)
+        approx(f; kw...) = approximate(f, domain[T], method(); kw...)
         f = x -> exp(3x);
         r = approx(f, tol=1e-5)
         @test !pass(f, r, pts, atol=1e-10)
@@ -65,7 +65,7 @@
     end
 
     @testset "Poles, zeros, residues" for T in (Float64,)
-        approx(f; kw...) = approximate(f, domain[T]; kw...)
+        approx(f; kw...) = approximate(f, domain[T], Barycentric(); kw...)
         f = z -> (z+1) * (z+2) / ((z+3) * (z+4))
         r = approx(f)
         pol = poles(r)
@@ -86,7 +86,7 @@
     end
 
     @testset "Vertical scaling in $T" for T in (Float64, Double64)
-        approx(f; kw...) = approximate(f, domain[T]; method=Barycentric, kw...)
+        approx(f; kw...) = approximate(f, domain[T], Barycentric(); kw...)
         pts = test_points[T]
         f = x -> T(10)^50*sinh(x); @test pass(f, approx(f), pts, rtol=2000*eps(T))
         f = x -> T(10)^(-50)*cosh(x); @test pass(f, approx(f), pts, rtol=2000*eps(T))
@@ -96,7 +96,7 @@
         T = Float64
         pts = test_points[T]
         tol = 2000*eps(T)
-        approx(f; kw...) = approximate(f, domain[T]; method=Barycentric, kw...)
+        approx(f; kw...) = approximate(f, domain[T], Barycentric(); kw...)
         f = x -> 0; @test pass(f, approx(f), pts, atol=2e-13)
         f = x -> x; @test pass(f, approx(f), pts, atol=2e-13)
         f = x -> 1im*x; @test pass(f, approx(f), pts, atol=2e-13)
diff --git a/test/operations.jl b/test/operations.jl
index b23cf648..50a9d37e 100644
--- a/test/operations.jl
+++ b/test/operations.jl
@@ -13,7 +13,7 @@
                 (x -> log(1.1 - x), x -> -1 / (1.1 - x), x -> -1 / (1.1 - x)^2),
                 (sin, cos, x -> -sin(x)),
             ))
-                r = approximate(f, domain; method)
+                r = approximate(f, domain, method())
                 @test isapprox(derivative(r; allowed=true), df, atol=sqrt(eps()))
                 @test isapprox(derivative(r, 2; allowed=true), d2f, atol=50sqrt(eps()))
                 vals = derivative(r.fun, 0:2)(0.25)
@@ -26,12 +26,12 @@
 
     @testset "Arithmetic with $method" verbose=true for method in (Barycentric, Thiele)
         @testset "Domain $iter" for (iter, domain) in enumerate((unit_interval, Shapes.square))
-            e = approximate(exp, domain; method)
-            t = approximate(tan, domain; method)
+            e = approximate(exp, domain, method())
+            t = approximate(tan, domain, method())
             @test (e / e) ≈ 1
             @test (3im * t - 2im * t) ≈ 1im * t
             c = cis
-            ec = approximate(exp, unit_circle; method)
+            ec = approximate(exp, unit_circle, method())
             @testset "$(op)" for op in (+, -, *, /)
                 @test values(op(e, 3.14im)) ≈ op.(values(e), 3.14im)
                 @test nodes(op(e, 3.14im)) ≈ nodes(e)
@@ -49,7 +49,7 @@
     end
 
     @testset "Arithmetic with zero for $method" verbose=true for method in (Barycentric, Thiele)
-        r = approximate(exp, unit_interval; method)
+        r = approximate(exp, unit_interval, method())
         @test r + 0 ≈ r
         @test r - 0 ≈ r
         @test r * 0 ≈ 0
diff --git a/test/real_interval.jl b/test/real_interval.jl
index b993aba2..936a833f 100644
--- a/test/real_interval.jl
+++ b/test/real_interval.jl
@@ -11,7 +11,7 @@
         T = Float64
         tol = 3000*eps(T)
         pts = test_points[T]
-        approx(f; kw...) = approximate(f, domain[T]; method, kw...)
+        approx(f; kw...) = approximate(f, domain[T], method(); kw...)
         @testset "Function $iter" for (iter, f) in enumerate((
             exp,
             cis,
@@ -41,10 +41,10 @@
         T = Float64
         pts = test_points[T]
         f = x -> sin(40x) * exp(-8x^2)
-        r = approximate(f, domain[T], method=Thiele)
+        r = approximate(f, domain[T], Thiele())
         @test isapprox(f.(pts), r(pts), norm=u->maximum(abs, u), rtol=2e-11)
         f = x -> cis(16x)
-        r = approximate(f, domain[T], method=Thiele)
+        r = approximate(f, domain[T], Thiele())
         @test isapprox(f.(pts), r(pts), norm=u->maximum(abs, u), rtol=2e-11)
     end
 
@@ -52,7 +52,7 @@
         T = Double64
         tol = 2000*eps(T)
         pts = test_points[T]
-        approx(f; kw...) = approximate(f, domain[T]; method, kw...)
+        approx(f; kw...) = approximate(f, domain[T], method(); kw...)
         @testset "Function $iter" for (iter, f) in enumerate((
             x -> cis(x),
             x -> exp(x),
@@ -70,12 +70,12 @@
 
     @testset "Float type conversion for $method" for method in (Barycentric, Thiele)
         f = z -> abs(z - 1im)
-        r = approximate(f, unit_interval; method)
+        r = approximate(f, unit_interval, method())
         r32 = convert(Float32, r.fun)
         @test r32 isa method{Float32,Float32}
 
         f = z -> cis(z)
-        r = approximate(f, unit_interval; method)
+        r = approximate(f, unit_interval, method())
         r32 = convert(Float32, r.fun)
         @test r32 isa method{Float32,ComplexF32}
     end
@@ -84,17 +84,60 @@
         T = Float64
         method = Barycentric
         pts = test_points[T]
-        approx(f; kw...) = approximate(f, domain[T]; method, kw...)
+        approx(f; kw...) = approximate(f, domain[T], method(); kw...)
         f = x -> exp(3x);
         r = approx(f, tol=1e-5)
         @test !pass(f, r, pts, atol=1e-7)
         @test pass(f, r, pts, atol=5e-5)
-        f = x -> abs(x);  @test pass(f, approx(f, stagnation=30), pts, atol=1e-10)
-        f = x -> abs(x - 0.95);  @test pass(f, approx(f, stagnation=30), pts, atol=1e-9)
+        # `abs` converges root-exponentially and needs degree ~116 to reach these
+        # tolerances, so it is given a budget larger than the default max_degree.
+        f = x -> abs(x);  @test pass(f, approx(f, stagnation=30, max_degree=150), pts, atol=1e-10)
+        f = x -> abs(x - 0.95);  @test pass(f, approx(f, stagnation=30, max_degree=150), pts, atol=1e-9)
+    end
+
+    @testset "Convergence status" begin
+        T = Float64
+        pts = test_points[T]
+
+        # An easy function reaches the tolerance.
+        for method in (Barycentric, Thiele)
+            r = approximate(exp, domain[T], method())
+            s = status(r)
+            @test s isa ConvergenceStatus
+            @test s.reason == :converged
+            @test isconverged(r) && isconverged(s)
+            @test s.iterations == length(get_history(r)[1])
+            @test s.error == get_history(r)[2][s.best]
+        end
+
+        # `abs` cannot reach the tolerance within the default degree budget, so the
+        # iteration is truncated rather than converged.
+        f = x -> abs(x)
+        r = approximate(f, domain[T], Barycentric(); stagnation=30, max_degree=40)
+        @test status(r).reason == :max_degree
+        @test !isconverged(r)
+        @test status(r).iterations == 40
+        # Raising the budget lets the same problem converge.
+        @test isconverged(approximate(f, domain[T], Barycentric(); stagnation=30, max_degree=150))
+
+        # A rewound approximation did not stop for the reason the iteration did.
+        r = approximate(exp, domain[T], Barycentric())
+        @test status(rewind(r, 3)).reason == :rewound
+        @test status(rewind(r, 3)).best == 3
+        @test !isconverged(rewind(r, 3))
+
+        # Discrete domains report status too.
+        r = approximate(exp, pts, Thiele())
+        @test isconverged(r)
+
+        # Prescribed poles are a direct solve, so there is no iteration to report.
+        r = approximate(x -> 1 / (x^2 + 4), domain[T], [2im, -2im])
+        @test isnothing(status(r))
+        @test !isconverged(r)
     end
 
     @testset "Nodes, values, degree for Barycentric" begin
-        r = approximate(exp, unit_interval; method=Barycentric)
+        r = approximate(exp, unit_interval, Barycentric())
         @test length(nodes(r)) == 6
         @test length(weights(r)) == 6
         @test minimum(nodes(r)) ≈ -1
@@ -106,14 +149,14 @@
         deg, err, zp, allowed, best = get_history(r)
         @test deg[end] == degree(r)
         @test length(deg) == length(err) == length(allowed)
-        r = approximate(exp, unit_interval; method=Barycentric, allowed=true)
+        r = approximate(exp, unit_interval, Barycentric(); allowed=true)
         deg, err, zp, allowed, best = get_history(r)
         @test deg[end] == degree(r)
         @test length(deg) == length(err) == length(allowed)
     end
 
     @testset "Nodes, values, degree for Thiele" begin
-        r = approximate(exp, unit_interval; method=Thiele)
+        r = approximate(exp, unit_interval, Thiele())
         @test length(nodes(r)) == 11
         @test length(weights(r)) == 11
         @test minimum(nodes(r)) ≈ -1
@@ -125,14 +168,14 @@
         deg, err, zp, allowed, best = get_history(r)
         @test deg[end] == degree(r)
         @test length(deg) == length(err) == length(allowed)
-        r = approximate(exp, unit_interval; method=Barycentric, allowed=true)
+        r = approximate(exp, unit_interval, Barycentric(); allowed=true)
         deg, err, zp, allowed, best = get_history(r)
         @test deg[end] == degree(r)
         @test length(deg) == length(err) == length(allowed)
     end
 
     @testset "Poles, zeros, residues in $T for Barycentric" for T in (Float64, Double64)
-        approx(f; kw...) = approximate(f, domain[T]; method=Barycentric, kw...)
+        approx(f; kw...) = approximate(f, domain[T], Barycentric(); kw...)
         f = z -> (z+1) * (z+2) / ((z+3) * (z+4))
         r = approx(f)
         pol = poles(r)
@@ -156,7 +199,7 @@
     end
 
     @testset "Poles, zeros, residues for Thiele" begin
-        approx(f; kw...) = approximate(f, domain[Float64]; method=Thiele, kw...)
+        approx(f; kw...) = approximate(f, domain[Float64], Thiele(); kw...)
         f = z -> (z+1) * (z+2) / ((z+3) * (z+4))
         r = approx(f)
         pol = poles(r)
@@ -186,7 +229,7 @@
 
     @testset "Vertical scaling in $T" for T in (Float64, Double64)
         pts = test_points[T]
-        approx(f; kw...) = approximate(f, domain[T]; method=Barycentric, kw...)
+        approx(f; kw...) = approximate(f, domain[T], Barycentric(); kw...)
         f = x -> T(10)^50 * sin(x); @test pass(f, approx(f), pts, rtol=2000*eps(T))
         f = x -> T(10)^(-50) * cos(x); @test pass(f, approx(f), pts, rtol=2000*eps(T))
         # TODO: Horizontal scaling is broken, because pole computation uses 1
@@ -199,7 +242,7 @@
         T = Float64
         tol = 2000*eps(T)
         pts = test_points[T]
-        approx(f; kw...) = approximate(f, domain[T]; method=Barycentric, kw...)
+        approx(f; kw...) = approximate(f, domain[T], Barycentric(); kw...)
         f = x -> 0; @test pass(f, approx(f, max_iter=1), pts, atol=tol)
         f = x -> x; @test pass(f, approx(f, max_iter=2), pts, atol=tol)
         f = x -> 1im*x; @test pass(f, approx(f, max_iter=2), pts, atol=tol)
@@ -213,9 +256,34 @@
 
     @testset "Interval [$a, $b]" for (a, b) in ((-2, 3), (0, 4), (-2e-4, 0), (-3e3, 5e6))
         pts = range(a, b, 1000)
-        approx(f; kw...) = approximate(f, Segment(a, b); method=Barycentric, kw...)
+        approx(f; kw...) = approximate(f, Segment(a, b), Barycentric(); kw...)
         toler(f) = 1000 * eps() * max(b - a, abs(f(a)), abs(f(b)))
         f = x -> 1 / sin(a + (b-a)*1.05im - x); @test pass(f, approx(f), pts, atol=toler(f))
         f = x -> exp(-10/(a + 1.1*(b-a) - x)); @test pass(f, approx(f), pts, atol=toler(f))
     end
+
+    @testset "ClosedInterval from IntervalSets" begin
+        pts = test_points[Float64]
+        f = x -> exp(x)
+        r = approximate(f, -1..1)
+        @test RFA.domain(r) == Segment(-1.0, 1.0)
+        @test pass(f, r, pts, rtol=3000*eps())
+
+        r_bary = approximate(f, -1..1, Barycentric())
+        @test RFA.domain(r_bary) == Segment(-1.0, 1.0)
+        @test pass(f, r_bary, pts, rtol=3000*eps())
+
+        g = x -> 1 / (x^2 + 4)
+        r_pf = approximate(g, -1..1, [2im, -2im])
+        @test RFA.domain(r_pf) == Segment(-1.0, 1.0)
+        @test pass(g, r_pf, pts, rtol=3000*eps())
+
+        r_ab = approximate(f, 0.5..2.5)
+        @test RFA.domain(r_ab) == Segment(0.5, 2.5)
+    end
+
+    @testset "Single-argument call for -1..1" begin
+        r = approximate(exp)
+        @test RFA.domain(r) == Segment(-1.0, 1.0)
+    end
 end