Add Traub method in Newton solver

Author: facelessuser

coloraide/algebra.py

@@ -511,13 +511,13 @@ def solve_newton(
     maxiter: int = 50,
     rtol: float = RTOL,
     atol: float = ATOL,
-    ostrowski: bool = False,
+    order: int = 2,
     bounds: tuple[float, float] | None = None,
 ) -> 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.
+    If order <= 3 and second derivative is given, Halley's method will be used as an additional step.
     Newton provides 2nd order convergence and Halley provides 3rd order convergence.
 
     ```
@@ -535,7 +535,16 @@ def solve_newton(
     Step4: halley = xn - yn / (1 - 0.5 * yn * f''(xn) / f'(xn))
     ```
 
-    If Ostrowski method is enabled, only one derivative is needed, but you can get 4th order convergence.
+    If order == 3 and 2nd derivative is not provided, we can use Traub's method which gives 3rd order convergence
+    without a second derivative.
+
+    ```
+    yn = xn - f(xn) / f'(xn)
+    traub = yn - f(yn) / f'(xn)
+    ```
+
+    If order >= 3, we can use Ostrowski method where only one derivative is needed, but you can get 4th order
+    convergence.
 
     ```
     yn = xn - f(xn) / f'(xn)
@@ -545,6 +554,9 @@ def solve_newton(
     Return result along with True if converged, False if did not converge, None if could not converge.
     """
 
+    if dx2 is not None and order < 3:
+        order = 3
+
     if bounds is not None:
         lo, hi = bounds
         bracketed = True
@@ -579,22 +591,26 @@ def solve_newton(
             return x0, None
 
         # Newton step
-        newton =  fx / d1
-
-        # If second derivative is provided, apply the Halley's method step: 3rd order convergence.
-        # The Newton step has been factored out of Halley's such that we can apply the rest to make
-        # it Halley's, and if we can't, or shouldn't apply it, it remains a Newton step.
-        if dx2 is not None and not ostrowski:
-            d2 = dx2(x0, *args) if args else dx2(x0)
-            denom = 1 - (0.5 * newton * d2) / d1
-            if abs(denom) >= ATOL:
-                newton /= denom
+        newton = fx / d1
+
+        if order == 3:
+            # Halley's method: 3rd order convergence.
+            if dx2 is not None:
+                d2 = dx2(x0, *args) if args else dx2(x0)
+                denom = 1 - (0.5 * newton * d2) / d1
+                if abs(denom) >= ATOL:
+                    newton /= denom
+
+            # Traub's method: 3rd order convergence
+            else:
+                fy = f0(x0 - newton, *args) if args else f0(x0 - newton)
+                newton -= fy / d1
 
         # If change is under our epsilon, we can consider the result converged.
         x0 -= newton
 
-        # Use Ostrowski method: 4th order convergence
-        if ostrowski:
+        # Ostrowski's method: 4th order convergence
+        if order == 4:
             fy = f0(x0, *args) if args else f0(x0)
             denom = fx - 2 * fy
             if abs(denom) >= ATOL:

docs/src/dictionary/en-custom.txt

@@ -232,6 +232,7 @@ TODO
 TORTIOUS
 Tetradic
 Tokenize
+Traub
 Triadic
 Trichromacy
 Trilinear

docs/src/markdown/colors/ipt.md

@@ -31,7 +31,7 @@ The IPT color appearance model excels at providing a formulation for hue where a
 perceived hue independent of the values of lightness and chroma (which is the general ideal for any color appearance
 model, but hard to achieve). It is therefore well-suited for gamut mapping implementations.
 
-[Learn more](https://www.researchgate.net/publication/21677980_Development_and_Testing_of_a_Color_Space_IPT_with_Improved_Hue_Uniformity.).
+[Learn more](https://www.researchgate.net/publication/21677980_Development_and_Testing_of_a_Color_Space_IPT_with_Improved_Hue_Uniformity).
 ///
 
 ## Channel Aliases

tests/test_algebra.py

@@ -2851,8 +2851,8 @@ def test_sign(self):
         self.assertEqual(alg.sign(float('inf')), 1)
         self.assertEqual(alg.sign(float('-inf')), -1)
 
-    def test_solve_newton_and_halley(self):
-        """Test solving with Newton."""
+    def test_solve_newton_variants(self):
+        """Test solving with Newton variants."""
 
         def f0(x):
             return 2 * x ** 3 - 3 * x ** 2 - 3 * x + 2
@@ -2869,8 +2869,10 @@ def dx2(x):
         self.assertEqual(alg.solve_newton(1, f0, dx), (0.5, True))
         # Halley
         self.assertEqual(alg.solve_newton(1, f0, dx, dx2), (0.5, True))
+        # Traub
+        self.assertEqual(alg.solve_newton(1, f0, dx, order=3), (0.5, True))
         # Ostrowski
-        self.assertEqual(alg.solve_newton(1, f0, dx, ostrowski=True), (0.5, True))
+        self.assertEqual(alg.solve_newton(1, f0, dx, order=4), (0.5, True))
 
 
     def test_solve_newton_bracketed(self):