python · PhysicistJohn · Jul 24, 2026 · Jul 24, 2026
diff --git a/Doc/library/math.rst b/Doc/library/math.rst
@@ -36,6 +36,7 @@ noted otherwise, all return values are floats.
 :func:`fmax(x, y) <fmax>`                             Maximum of two floating-point values
 :func:`fmin(x, y) <fmin>`                             Minimum of two floating-point values
 :func:`fmod(x, y) <fmod>`                             Remainder of division ``x / y``
+:func:`lerp(a, b, t) <lerp>`                          Linear interpolation between *a* and *b* using parameter *t*
 :func:`modf(x) <modf>`                                Fractional and integer parts of *x*
 :func:`remainder(x, y) <remainder>`                   Remainder of *x* with respect to *y*
 :func:`trunc(x) <trunc>`                              Integer part of *x*
@@ -351,6 +352,31 @@ Floating point manipulation functions
    :func:`frexp`.
 
 
+.. function:: lerp(a, b, t)
+
+   Linear interpolation between *a* and *b*, using *t* as the interpolation
+   parameter. Returns ``a + t * (b - a)``, computed in a way that is
+   guaranteed to be numerically well-behaved:
+
+   * ``lerp(a, b, 0.0)`` returns *a* exactly.
+   * ``lerp(a, b, 1.0)`` returns *b* exactly.
+   * ``lerp(a, b, t)`` is monotonic in *t*, for finite *a* and *b* with
+     ``a <= b`` and finite *t*.
+   * *t* is not restricted to ``[0.0, 1.0]``; values outside that range
+     perform extrapolation.
+
+   The naive formula ``a + t * (b - a)`` does not guarantee any of the
+   above; it can, for example, fail to return *b* exactly when *t* is
+   ``1.0`` due to floating-point rounding. :func:`lerp` uses the
+   algorithm standardized for C++20's ``std::lerp``, which fixes these
+   issues. As a consequence, results involving infinities or NaN follow
+   from ordinary floating-point arithmetic rules on each branch of that
+   algorithm, rather than a general guarantee that NaN inputs always
+   produce a NaN result.
+
+   .. versionadded:: 3.16
+
+
 .. function:: nextafter(x, y, steps=1)
 
    Return the floating-point value *steps* steps after *x* towards *y*.

@@ -371,6 +371,11 @@ math
   754-2019 and standardized in C23.
   (Contributed by Jeff Epler in :gh:`150534`.)
 
+* Added :func:`math.lerp` for linear interpolation, using the algorithm
+  standardized for C++20's ``std::lerp``, which guarantees exact
+  endpoints and monotonicity.
+  (Contributed by PhysicistJohn in :gh:`154573`.)
+
 
 os
 --

diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py
@@ -2847,6 +2847,112 @@ def assertIsNegativeZero(self, value):
         )
 
 
+class LerpTests(unittest.TestCase):
+    """ Tests for math.lerp. """
+
+    def test_basic(self):
+        self.assertEqual(math.lerp(0.0, 10.0, 0.5), 5.0)
+        self.assertEqual(math.lerp(1.0, 3.0, 2.0), 5.0)  # extrapolation
+        self.assertEqual(math.lerp(1.0, 3.0, -1.0), -1.0)  # extrapolation
+        self.assertEqual(math.lerp(-5.0, 5.0, 0.5), 0.0)
+
+    def test_int_arguments_are_accepted(self):
+        # Clinic converts ints to float, like other math functions.
+        self.assertEqual(math.lerp(0, 10, 0), 0.0)
+        self.assertEqual(math.lerp(0, 10, 1), 10.0)
+
+    def test_endpoints_are_exact(self):
+        # t == 0.0 and t == 1.0 must return the endpoints exactly, even
+        # where the naive "a + t * (b - a)" formula would round off.
+        cases = [
+            (0.0, 10.0), (10.0, 0.0), (-5.0, 5.0), (1e300, -1e300),
+            (524560164915884.0, -995787893297778.6),
+        ]
+        for a, b in cases:
+            with self.subTest(a=a, b=b):
+                self.assertEqual(math.lerp(a, b, 0.0), a)
+                self.assertEqual(math.lerp(a, b, 1.0), b)
+
+    def test_naive_formula_would_be_inexact_at_t1(self):
+        # Concrete demonstration of why lerp() special-cases t == 1.0
+        # rather than using "a + t * (b - a)" directly.
+        a = 524560164915884.0
+        b = -995787893297778.6
+        naive = a + 1.0 * (b - a)
+        self.assertNotEqual(naive, b)
+        self.assertEqual(math.lerp(a, b, 1.0), b)
+
+    def test_a_equals_b(self):
+        # lerp(a, a, t) == a for any finite t, including outside [0, 1].
+        for t in [-100.0, -1.0, 0.0, 0.5, 1.0, 2.0, 100.0]:
+            with self.subTest(t=t):
+                self.assertEqual(math.lerp(3.5, 3.5, t), 3.5)
+
+    def test_monotonic(self):
+        random.seed(163875206)
+        for _ in range(1000):
+            a = random.uniform(-1e10, 1e10)
+            b = random.uniform(-1e10, 1e10)
+            lo, hi = min(a, b), max(a, b)
+            ts = sorted(random.uniform(-2.0, 2.0) for _ in range(8))
+            values = [math.lerp(lo, hi, t) for t in ts]
+            self.assertEqual(values, sorted(values))
+
+    def test_matches_cpp20_stdlerp_reference_values(self):
+        # Bit-for-bit reference values captured from libc++'s std::lerp
+        # (LLVM main, clang -std=c++20), including edge cases involving
+        # zero, signed infinities, and NaN. math.lerp intentionally
+        # implements the same standardized algorithm, so it should match
+        # exactly, including cases where NaN in 'a' does not propagate
+        # (a documented characteristic of the algorithm itself, not a
+        # special case added here).
+        nan = math.nan
+        inf = math.inf
+        reference = [
+            ((0.0, 10.0, 0.5), 5.0),
+            ((10.0, 0.0, 0.0), 10.0),
+            ((10.0, 0.0, 1.0), 0.0),
+            ((1.0, 1.0, 100.0), 1.0),
+            ((1.0, 1.0, -100.0), 1.0),
+            ((-0.0, 10.0, 0.5), 5.0),
+            ((1e300, 1e300, 5.0), 1e300),
+            ((1e308, -1e308, 3.0), -inf),
+            ((-1e308, 1e308, 3.0), inf),
+            ((5.0, -5.0, 0.5), 0.0),
+        ]
+        for (a, b, t), expected in reference:
+            with self.subTest(a=a, b=b, t=t):
+                self.assertEqual(math.lerp(a, b, t), expected)
+
+        nan_cases = [
+            (5.0, nan, 0.5),
+            (0.0, 10.0, nan),
+            (-inf, inf, 0.5),
+        ]
+        for a, b, t in nan_cases:
+            with self.subTest(a=a, b=b, t=t):
+                self.assertTrue(math.isnan(math.lerp(a, b, t)))
+
+        # These do NOT produce NaN, even though 'a' or 't' is infinite/NaN
+        # in one operand: this matches std::lerp's documented behavior,
+        # where the t == 1 and clamping branches can return a or b
+        # without evaluating an expression that would propagate NaN.
+        self.assertEqual(math.lerp(nan, 5.0, 0.5), 5.0)
+        self.assertEqual(math.lerp(inf, 5.0, 0.5), 5.0)
+        self.assertEqual(math.lerp(5.0, inf, 0.5), inf)
+
+    def test_wrong_number_of_arguments(self):
+        self.assertRaises(TypeError, math.lerp)
+        self.assertRaises(TypeError, math.lerp, 1.0)
+        self.assertRaises(TypeError, math.lerp, 1.0, 2.0)
+        self.assertRaises(TypeError, math.lerp, 1.0, 2.0, 3.0, 4.0)
+
+    def test_wrong_argument_type(self):
+        self.assertRaises(TypeError, math.lerp, 'a', 2.0, 0.5)
+        self.assertRaises(TypeError, math.lerp, 1.0, 'b', 0.5)
+        self.assertRaises(TypeError, math.lerp, 1.0, 2.0, 't')
+
+
 def load_tests(loader, tests, pattern):
     from doctest import DocFileSuite
     tests.addTest(DocFileSuite(os.path.join("mathdata", "ieee754.txt")))

diff --git a/Misc/NEWS.d/next/Library/2026-07-23-23-00-00.gh-issue-154573.math-lerp.rst b/Misc/NEWS.d/next/Library/2026-07-23-23-00-00.gh-issue-154573.math-lerp.rst
@@ -0,0 +1,5 @@
+Add :func:`math.lerp`, computing linear interpolation between two
+values using the algorithm standardized for C++20's ``std::lerp``,
+which guarantees exact endpoints (``lerp(a, b, 0.0) == a`` and
+``lerp(a, b, 1.0) == b``) and monotonicity, properties the naive
+``a + t * (b - a)`` formula does not have. Patch by PhysicistJohn.
diff --git a/Modules/clinic/mathmodule.c.h b/Modules/clinic/mathmodule.c.h
diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c
@@ -1719,6 +1719,53 @@ math_ldexp_impl(PyObject *module, double x, PyObject *i)
 }
 
 
+/*[clinic input]
+math.lerp
+
+    a: double
+    b: double
+    t: double
+    /
+
+Linear interpolation between a and b using parameter t.
+
+Returns a + t * (b - a), computed in a way that guarantees intended
+mathematical behavior:
+
+* lerp(a, b, 0.0) == a
+* lerp(a, b, 1.0) == b
+* lerp(a, b, t) is monotonic in t, when a <= b
+* t outside [0.0, 1.0] is allowed, for extrapolation
+[clinic start generated code]*/
+
+static PyObject *
+math_lerp_impl(PyObject *module, double a, double b, double t)
+/*[clinic end generated code: output=7ca02ba2e81fd287 input=978d065447b6b7e0]*/
+{
+    /* This is the algorithm standardized for std::lerp in C++20
+     * ([c.math.lerp], from WG21 P0811R3 "Well-behaved interpolation
+     * for numbers and pointers" by S. Davis Herring). It guarantees
+     * properties the naive "a + t * (b - a)" formula does not: exact
+     * endpoints, monotonicity, and lerp(a, a, t) == a for finite t.
+     */
+    if ((a <= 0.0 && b >= 0.0) || (a >= 0.0 && b <= 0.0)) {
+        return PyFloat_FromDouble(t * b + (1.0 - t) * a);
+    }
+
+    if (t == 1.0) {
+        return PyFloat_FromDouble(b);
+    }
+
+    double x = a + t * (b - a);
+    if ((t > 1.0) == (b > a)) {
+        return PyFloat_FromDouble(b < x ? x : b);
+    }
+    else {
+        return PyFloat_FromDouble(x < b ? x : b);
+    }
+}
+
+
 /*[clinic input]
 math.modf
 
@@ -3264,6 +3311,7 @@ static PyMethodDef math_methods[] = {
     MATH_ISINF_METHODDEF
     MATH_ISNAN_METHODDEF
     MATH_LDEXP_METHODDEF
+    MATH_LERP_METHODDEF
     {"lgamma",          math_lgamma,    METH_O,         math_lgamma_doc},
     {"log",             _PyCFunction_CAST(math_log),       METH_FASTCALL,  math_log_doc},
     {"log1p",           math_log1p,     METH_O,         math_log1p_doc},