Inverse quadratic interpolation in Python
Fit a quadratic in x as a function of y through the last three points and evaluate it at y = 0. One evaluation in each step and order 1.839 at a simple root, but no bracket.
- File
inverse_quadratic_interpolation.py, 69 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/inverse-quadratic-interpolation/. The authors are listed in prov/6ZVEEE.jsonl.
"""Inverse quadratic interpolation: fit x as a quadratic in y through three points, take y = 0."""
import math
from dataclasses import dataclass, field
from typing import Callable, List
@dataclass
class IQIResult:
root: float
iterations: int
evaluations: int
status: str # "step", "zero", "equal-values" or "max-iter"
iterates: List[float] = field(default_factory=list)
@property
def converged(self) -> bool:
return self.status in ("step", "zero")
def iqi_step(x0: float, f0: float, x1: float, f1: float, x2: float, f2: float) -> float:
"""The zero of the quadratic in y through (f0, x0), (f1, x1), (f2, x2), in Newton form.
Precondition: f0, f1, f2 are distinct.
"""
d21 = (x2 - x1) / (f2 - f1)
d10 = (x1 - x0) / (f1 - f0)
d210 = (d10 - d21) / (f0 - f2)
return x2 - f2 * d21 + f1 * f2 * d210
def inverse_quadratic_interpolation(f: Callable[[float], float], x0: float, x1: float, x2: float,
tol: float = 1e-12, max_iter: int = 100) -> IQIResult:
"""Iterate inverse quadratic interpolation from three starting points.
Precondition: x0, x1, x2 near a simple root of f, with distinct values
of f. Local method: no bracket and no guarantee from far away. Stops when
a step is at most tol (a heuristic test), when f is exactly zero, when two
of the three values of f are equal, or after max_iter steps. Cost: one
evaluation of f in each step plus three at the start; order about 1.839
near a simple root. Raises ArithmeticError if f returns NaN or inf.
"""
f0, f1, f2 = f(x0), f(x1), f(x2) # @line L1
evals = 3
iterates: List[float] = []
for v in (f0, f1, f2):
if math.isnan(v) or math.isinf(v):
raise ArithmeticError("f is not finite at a starting point")
for x, v in ((x0, f0), (x1, f1), (x2, f2)):
if v == 0:
return IQIResult(x, 0, evals, "zero", iterates)
for k in range(1, max_iter + 1): # @line L2
if f0 == f1 or f1 == f2 or f0 == f2: # @line L3
return IQIResult(x2, k - 1, evals, "equal-values", iterates)
x3 = iqi_step(x0, f0, x1, f1, x2, f2) # @line L4
f3 = f(x3) # @line L5
evals += 1
if math.isnan(f3) or math.isinf(f3):
raise ArithmeticError("f is not finite at an iterate")
iterates.append(x3)
if f3 == 0: # @line L6
return IQIResult(x3, k, evals, "zero", iterates)
if abs(x3 - x2) <= tol:
return IQIResult(x3, k, evals, "step", iterates)
x0, f0, x1, f1, x2, f2 = x1, f1, x2, f2, x3, f3 # @line L7
return IQIResult(x2, max_iter, evals, "max-iter", iterates)
test_inverse_quadratic_interpolation.py 121 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/inverse-quadratic-interpolation/. The authors are listed in prov/6ZVEEE.jsonl.
import math
import random
import unittest
from decimal import Decimal, getcontext
from inverse_quadratic_interpolation import inverse_quadratic_interpolation, iqi_step
CUBIC_ROOT = 1.5213797068045676
Q = 1.839286755214161 # real root of q^3 = q^2 + q + 1
class Counter:
def __init__(self, f):
self.f, self.n = f, 0
def __call__(self, x):
self.n += 1
return self.f(x)
class TestIQI(unittest.TestCase):
def test_cubic_table(self):
res = inverse_quadratic_interpolation(lambda x: x ** 3 - x - 2, 1.0, 1.5, 2.0, 1e-14)
self.assertTrue(res.converged)
self.assertAlmostEqual(res.iterates[0], 1.5272727272727273, places=14)
self.assertAlmostEqual(res.iterates[1], 1.521336540179002, places=12)
self.assertLess(abs(res.root - CUBIC_ROOT), 1e-14)
def test_step_is_exact_on_inverse_quadratic(self):
# f is the inverse of x = 2 + 3y + y^2 near y = 0, so one step is exact
g = lambda y: 2 + 3 * y + y * y
ys = [0.1, -0.2, 0.3]
x = iqi_step(g(ys[0]), ys[0], g(ys[1]), ys[1], g(ys[2]), ys[2])
self.assertAlmostEqual(x, 2.0, places=14)
def test_evaluation_count(self):
f = Counter(lambda x: math.cos(x) - x)
res = inverse_quadratic_interpolation(f, 0.0, 0.5, 1.0, 1e-12)
self.assertEqual(f.n, res.evaluations)
self.assertEqual(res.evaluations, res.iterations + 3)
self.assertLess(abs(res.root - 0.7390851332151607), 1e-15)
def test_error_constant_and_order(self):
getcontext().prec = 120
f = lambda x: x ** 3 - x - 2
r = Decimal("1.5")
for _ in range(20):
r = r - f(r) / (3 * r * r - 1)
fp, fpp = 3 * r * r - 1, 6 * r
C = (3 * fpp ** 2 - fp * 6) / (6 * fp ** 2)
xs = [Decimal(1), Decimal("1.5"), Decimal(2)]
fs = [f(x) for x in xs]
for _ in range(6):
x0, x1, x2 = xs[-3:]
f0, f1, f2 = fs[-3:]
x3 = iqi_step(x0, f0, x1, f1, x2, f2)
xs.append(x3)
fs.append(f(x3))
es = [abs(x - r) for x in xs]
self.assertAlmostEqual(float(es[8] / (es[7] * es[6] * es[5])), float(C), places=3)
logs = [float(es[i].ln() / es[i - 1].ln()) for i in range(6, 9)]
self.assertAlmostEqual(sum(logs) / len(logs), Q, delta=0.08)
def test_known_functions(self):
for f, s, r in ((lambda x: math.exp(x) - 2, (0.0, 0.5, 1.0), math.log(2)),
(lambda x: (x - 0.25) * (x - 3) * (x + 4), (0.0, 0.2, 0.5), 0.25)):
res = inverse_quadratic_interpolation(f, *s, 1e-13)
self.assertTrue(res.converged)
self.assertLess(abs(res.root - r), 1e-13)
self.assertLess(res.evaluations, 15)
def test_double_root_is_slow(self):
res = inverse_quadratic_interpolation(lambda x: (x - 1) ** 2 * (x + 2), 0.0, 0.5, 0.8, 1e-8, 200)
self.assertLess(abs(res.root - 1.0), 1e-4)
self.assertGreater(res.iterations, 15)
def test_divergence_far_from_root(self):
res = inverse_quadratic_interpolation(math.atan, 2.0, 3.0, 4.0, 1e-12, 6)
self.assertFalse(res.converged)
self.assertGreater(abs(res.iterates[-1]), 1e9)
def test_equal_values(self):
res = inverse_quadratic_interpolation(lambda x: x * x - 1, -2.0, 2.0, 0.5)
self.assertEqual(res.status, "equal-values")
self.assertEqual(res.evaluations, 3)
def test_zero_at_start_and_iterate(self):
res = inverse_quadratic_interpolation(lambda x: x - 1.0, 1.0, 2.0, 3.0)
self.assertEqual((res.root, res.status), (1.0, "zero"))
res = inverse_quadratic_interpolation(lambda x: 2 * x - 1, 0.0, 0.25, 1.0)
self.assertEqual((res.root, res.status), (0.5, "zero"))
def test_nan(self):
with self.assertRaises(ArithmeticError):
inverse_quadratic_interpolation(lambda x: 1.0 if x < 0.5 else math.nan, 0.0, 0.2, 0.7)
with self.assertRaises(ArithmeticError):
inverse_quadratic_interpolation(lambda x: math.inf if x > 2.5 else x - 3, 0.0, 1.0, 2.0)
def test_max_iter(self):
res = inverse_quadratic_interpolation(lambda x: x ** 3 - x - 2, 1.0, 1.5, 2.0, 0.0, 2)
self.assertEqual(res.status, "max-iter")
self.assertEqual(res.evaluations, 5)
def test_random(self):
rng = random.Random(1839)
for _ in range(600):
r = rng.uniform(-3, 3)
o1, o2 = r + rng.uniform(3, 6), r - rng.uniform(3, 6)
f = lambda x, r=r, o1=o1, o2=o2: (x - r) * (x - o1) * (x - o2)
h = rng.uniform(0.05, 0.3)
res = inverse_quadratic_interpolation(f, r - h, r + h / 3, r + h, 1e-12)
self.assertTrue(res.converged, msg=str((r, o1, o2, h, res.status)))
self.assertLess(abs(res.root - r), 1e-10)
self.assertLess(res.evaluations, 20)
if __name__ == "__main__":
unittest.main()
What the next level needs
- the tests did not run in this build