Better approximation of Bezier with Newton

Fix up our Newton approach
- Generalize Newton approach (optionally extended for Halley's method)
- Because we are setting it equal to our target, we need to check for
  zero.
- Kick out when our change is below the epsilon

coloraide/algebra.py

@@ -209,6 +209,66 @@ def polar_to_rect(c: float, h: float) -> tuple[float, float]:
     return a, b
 
 
+def solve_newton(
+    x0: float,
+    f0: Callable[[float], float],
+    dx: Callable[[float], float],
+    dx2: Callable[[float], float] | None = None,
+    maxiter: int = 8,
+    epsilon: float = 1e-12
+) -> tuple[float, bool | None]:
+    """
+    Solve equation using Newton's method.
+
+    If the second derivative is given, Halley's method will be used as an additional step.
+
+    ```
+    newton = f0 / dx
+    halley = (f0 * dx) / (dx ** 2 - 0.5 * f0 * dx)
+    ```
+
+    Algebraically, we can pull the Newton stop out of the Halley method into two separate steps
+    that can be applied on top of each other.
+
+    ```
+    Step1: newton = f0 / dx
+    Step2: halley = newton * (1 - 0.5 * newton * d2 / d1)
+    ```
+    """
+
+    for _ in range(maxiter):
+        # Get result form equation when setting value to expected result
+        f = f0(x0)
+
+        # If the result is zero, we've converged
+        if f == 0:
+            return x0, True
+
+        # Cannot find a solution if derivative is zero
+        d1 = dx(x0)
+        if d1 == 0:
+            return x0, None
+
+        # Calculate new, hopefully closer value with Newton's method
+        newton =  f / d1
+
+        # If second derivative is provided, apply the additional Halley's method step
+        if dx2 is not None:
+            d2 = dx2(x0)
+            value = (0.5 * newton * d2) / d1
+            # If the value is greater than one, the solution is deviating away from the newton step
+            if abs(value) < 1:
+                newton *= 1 - value
+
+        # If change is under our epsilon, we can consider the result converged.
+        prev = x0
+        x0 -= newton
+        if abs(x0 - prev) < epsilon:
+            return x0, True
+
+    return x0, False
+
+
 ################################
 # Interpolation and splines
 ################################

coloraide/easing.py

@@ -40,28 +40,28 @@
 """
 from __future__ import annotations
 import functools
-import math
+from . import algebra as alg
 from typing import Callable
 
 EPSILON = 1e-6
 MAX_ITER = 8
 
 
-def _bezier(t: float, a: float, b: float, c: float) -> float:
+def _bezier(a: float, b: float, c: float, y: float = 0.0) -> Callable[[float], float]:
     """
     Calculate the bezier point.
 
     We know that P0 and P3 are always (0, 0) and (1, 1) respectively.
     Knowing this we can simplify the equation by precalculating them in.
     """
 
-    return a * t ** 3 + b * t ** 2 + c * t
+    return lambda x: a * x ** 3 + b * x ** 2 + c * x - y
 
 
-def _bezier_derivative(t: float, a: float, b: float, c: float) -> float:
+def _bezier_derivative(a: float, b: float, c: float) -> Callable[[float], float]:
     """Derivative of curve."""
 
-    return 3 * a * t ** 2 + 2 * b * t + c
+    return lambda x: 3 * a * x ** 2 + 2 * b * x + c
 
 
 def _solve_bezier(target: float, a: float, b: float, c: float) -> float:
@@ -75,41 +75,27 @@ def _solve_bezier(target: float, a: float, b: float, c: float) -> float:
     # Try Newtons method to see if we can find a suitable value
     x = 0.0
     t = 0.5
-    last = math.nan
-    for _ in range(MAX_ITER):
-        # See how close we are to the desired `x`
-        x = _bezier(t, a, b, c) - target
-
-        # Get the derivative, but bail if it is too small,
-        # we will just keep looping otherwise.
-        dx = _bezier_derivative(t, a, b, c)
-
-        # Derivative is zero, we can't continue
-        if dx == 0:  # pragma: no cover
-            break
-
-        # Calculate new time and try again
-        t -= (x / dx)
-
-        # We converged on an solution
-        if t == last:  # pragma: no cover
-            return t
-
-        last = t
-
-    # We didn't fully converge but we are close, closer than our bisect epsilon
-    if abs(_bezier(t, a, b, c) - target) < EPSILON:
+    f0 = _bezier(a, b, c, y=target)
+    t, converged = alg.solve_newton(
+        t,
+        f0,
+        _bezier_derivative(a, b, c),
+        maxiter=MAX_ITER
+    )
+
+    # We converged or we are close enough
+    if converged or abs(t - target) < EPSILON:
         return t
 
     # We couldn't achieve our desired accuracy with Newton,
     # so bisect at lower accuracy until we arrive at a suitable value
     low, high = 0.0, 1.0
     t = target
-    while abs(high - low) > EPSILON:
-        x = _bezier(t, a, b, c)
-        if abs(x - target) < EPSILON:
+    while abs(high - low) >= EPSILON:
+        x = f0(t)
+        if abs(x) < EPSILON:
             return t
-        if x > target:
+        if x > 0:
             high = t
         else:
             low = t
@@ -187,7 +173,7 @@ def _calc_bezier(
     t = _solve_bezier(target, a[0], b[0], c[0])
 
     # Use the found `t` to locate the `y`
-    y = _bezier(t, a[1], b[1], c[1])
+    y = _bezier(a[1], b[1], c[1])(t)
     return y
 
 

coloraide/spaces/ryb.py

@@ -32,14 +32,14 @@ def srgb_to_ryb(rgb: Vector, cube_t: Matrix, cube: Matrix, biased: bool) -> Vect
     # Calculate the RYB value
     ryb = alg.ilerp3d(cube_t, rgb, vertices_t=cube)
     # Remove smoothstep easing if "biased" is enabled.
-    return [_solve_bezier(t, -2, 3, 0) if 0 <= t <= 1 else t for t in ryb] if biased else ryb
+    return [_solve_bezier(t, -2.0, 3.0, 0.0) if 0 <= t <= 1 else t for t in ryb] if biased else ryb
 
 
 def ryb_to_srgb(ryb: Vector, cube_t: Matrix, biased: bool) -> Vector:
     """Convert RYB to sRGB."""
 
     # Bias interpolation towards corners if "biased" enable. Bias is a smoothstep easing function.
-    return alg.lerp3d(cube_t, [_bezier(t, -2, 3, 0) if 0 <= t <= 1 else t for t in ryb] if biased else ryb)
+    return alg.lerp3d(cube_t, [_bezier(-2.0, 3.0, 0.0)(t) if 0 <= t <= 1 else t for t in ryb] if biased else ryb)
 
 
 class RYB(Regular, Space):

docs/src/markdown/about/changelog.md

@@ -1,5 +1,9 @@
 # Changelog
 
+## 4.2.2
+
+-   **FIX**: Speed up solving of cubic bezier for easing functions.
+
 ## 4.2.1
 
 -   **FIX**: Hex output should force gamut mapping even if it is requested to disable it as output will be invalid