simplify Laguerre solver

src/lib/Math/Bezier/Cubic.hs

@@ -182,7 +182,7 @@ closestPoint
 closestPoint pts c = pickClosest ( 0 :| 1 : roots ) -- todo: also include the self-intersection point if one exists
   where
     roots :: [ r ]
-    roots = filter ( \ r -> r > 0 && r < 1 ) ( realRoots 500 $ ddist @v pts c )
+    roots = filter ( \ r -> r > 0 && r < 1 ) ( realRoots 50 $ ddist @v pts c )
 
     pickClosest :: NonEmpty r -> ArgMin r ( r, p )
     pickClosest ( s :| ss ) = go s q nm0 ss

src/lib/Math/Bezier/Quadratic.hs

@@ -150,7 +150,7 @@ closestPoint
 closestPoint pts c = pickClosest ( 0 :| 1 : roots )
   where
     roots :: [ r ]
-    roots = filter ( \ r -> r > 0 && r < 1 ) ( realRoots 2000 $ ddist @v pts c ) 
+    roots = filter ( \ r -> r > 0 && r < 1 ) ( realRoots 50 $ ddist @v pts c ) 
 
     pickClosest :: NonEmpty r -> ArgMin r ( r, p )
     pickClosest ( s :| ss ) = go s q nm0 ss

src/lib/Math/Roots.hs

@@ -85,7 +85,8 @@ realRoots maxIters coeffs = mapMaybe isReal ( roots epsilon maxIters coeffs )
 -- Polynomial coefficients are given in order of descending degree (e.g. constant coefficient last).
 --
 -- N.B. The forward deflation process is only guaranteed to be numerically stable
--- if Laguerre's method finds roots in increasing order of magnitude.
+-- if Laguerre's method finds roots in increasing order of magnitude
+-- (which this function attempts to do).
 roots :: forall a. ( RealFloat a, Prim a, NFData a ) => a -> Int -> [ a ] -> [ Complex a ]
 roots eps maxIters coeffs = runST do
   let
@@ -97,10 +98,9 @@ roots eps maxIters coeffs = runST do
   let
     go :: Int -> [ Complex a ] -> ST s [ Complex a ]
     go !i rs = do
-      -- Estimate the root with minimum absolute value in order to
-      -- improve numerical stability of Laguerre's method with forward deflation.
-      !z0 <- minAbsStartPoint p
-      !r  <- force <$> laguerre eps maxIters p z0
+      -- Start at 0, attempting to find the root with smallest absolute value.
+      -- This improves numerical stability of the forward deflation process.
+      !r <- force <$> laguerre eps maxIters p 0
       if i <= 2
       then pure ( r : rs )
       else
@@ -114,41 +114,9 @@ roots eps maxIters coeffs = runST do
           go ( i - 2 ) ( r : conjugate r : rs )
   go sz []
 
--- | Estimate the root with smallest absolute value.
---
--- Polynomial coefficients are given in order of descending degree (e.g. constant coefficient last).
-minAbsStartPoint :: forall a m s. ( RealFloat a, Prim a, PrimMonad m, s ~ PrimState m ) => MutablePrimArray s a -> m ( Complex a )
-minAbsStartPoint p = do
-  n <- subtract 1 <$> getSizeofMutablePrimArray p
-  an <- readPrimArray p n
-  if abs an < epsilon
-  then pure 0
-  else do
-    a0 <- readPrimArray p 0
-    let
-      r, m0 :: Complex a
-      r  = log ( abs an :+ 0 )
-      m0 = exp $ ( r - log ( abs a0 :+ 0 ) ) / fromIntegral n
-      go :: Int -> Complex a -> m ( Complex a )
-      go i m
-        | i >= n
-        = pure m
-        | otherwise
-        = do
-          ai <- readPrimArray p i
-          if abs ai < epsilon
-          then go ( i + 1 ) m
-          else do
-            let
-              mi :: Complex a
-              mi = exp $ ( r - log ( abs ai :+ 0 ) ) / fromIntegral ( n - i )
-            if magnitude mi < magnitude m
-            then go ( i + 1 ) mi
-            else go ( i + 1 ) m
-    m <- go 1 m0
-    pure ( 0.5 * m )
-
 -- | Forward deflation of a polynomial by a root: factors out the root.
+--
+-- Mutates the input polynomial.
 deflate :: forall a m s. ( Num a, Prim a, PrimMonad m, s ~ PrimState m ) => a -> MutablePrimArray s a -> m ()
 deflate r p = do
   deg <- subtract 1 <$> getSizeofMutablePrimArray p
@@ -167,6 +135,8 @@ deflate r p = do
     go a0 1
 
 -- | Forward deflation of a polynomial with real coefficients by a pair of complex-conjugate roots.
+--
+-- Mutates the input polynomial.
 deflateConjugatePair :: forall a m s. ( Num a, Prim a, PrimMonad m, s ~ PrimState m ) => Complex a -> MutablePrimArray s a -> m ()
 deflateConjugatePair ( x :+ y ) p = do
   deg <- subtract 1 <$> getSizeofMutablePrimArray p
@@ -194,6 +164,8 @@ deflateConjugatePair ( x :+ y ) p = do
     go b1 a0 2
 
 -- | Laguerre's method.
+--
+-- Does not perform any mutation.
 laguerre
   :: forall a m s
   .  ( RealFloat a, Prim a, PrimMonad m, s ~ PrimState m )
@@ -215,6 +187,8 @@ laguerre eps maxIters p x0 = do
   go maxIters x0
 
 -- | Take a single step in Laguerre's method.
+--
+-- Does not perform any mutation.
 laguerreStep
   :: forall a m s
   .  ( RealFloat a, Prim a, PrimMonad m, s ~ PrimState m )
@@ -236,22 +210,33 @@ laguerreStep eps p p' p'' x = do
       g     = p'x / px
       g²    = g * g
       h     = g² - p''x / px
-      mult  = 1 --max 1 . min n . fromIntegral . round @_ @Integer . realPart $ log p'x / log ( px / p'x )
-      delta = sqrt $ ( n - 1 ) *: ( ( n / mult ) *: h - g² )
+    --m     = min n . max 1 . fromIntegral . sanitise . realPart $ negate ( log p'x / log g )
+      delta = sqrt $ ( n - 1 ) *: ( n *: h - g² ) -- sqrt $ ( ( n - m ) / m ) *: ( n *: h - g² )
       gp    = g + delta
       gm    = g - delta
       denom
         | magnitude gm > magnitude gp
-        = gm
+        = recip gm
         | otherwise
-        = gp
-    pure $ x - n *: ( recip denom )
+        = recip gp
+    pure $ x - n *: denom
 
   where
     (*:) :: a -> Complex a -> Complex a
     r *: (u :+ v) = ( r * u ) :+ ( r * v )
 
+{-
+    sanitise :: a -> Int
+    sanitise u
+      | isNaN u || isInfinite u
+      = 1
+      | otherwise
+      = round u
+-}
+
 -- | Evaluate a polynomial with real coefficients at a complex number.
+--
+-- Does not perform any mutation.
 eval
   :: forall a m s
   . ( RealFloat a, Prim a, PrimMonad m, s ~ PrimState m )
@@ -271,6 +256,8 @@ eval p x = do
   go ( an :+ 0 ) 1
 
 -- | Derivative of a polynomial.
+--
+-- Does not mutate its argument.
 derivative
   :: forall a m s
   . ( Num a, Prim a, PrimMonad m, s ~ PrimState m )