Modularise root isolation algorithms

Different root isolation algorithms now live in separate modules,
and are all instances of the RootIsolationAlgorithm typeclass.

This separates the algorithmic code from the top-level driver code
in Math.Root.Isolation.

MetaBrush.cabal

@@ -81,8 +81,9 @@ common common
         StandaloneDeriving
         StandaloneKindSignatures
         TupleSections
+        TypeAbstractions
         TypeApplications
-        TypeFamilies
+        TypeFamilyDependencies
         TypeOperators
         UnboxedTuples
         ViewPatterns

brush-strokes/brush-strokes.cabal

@@ -52,8 +52,9 @@ common common
         StandaloneDeriving
         StandaloneKindSignatures
         TupleSections
+        TypeAbstractions
         TypeApplications
-        TypeFamilies
+        TypeFamilyDependencies
         TypeOperators
         UnboxedTuples
         ViewPatterns
@@ -125,6 +126,10 @@ library
       , Math.Ring
       , Math.Roots
       , Math.Root.Isolation
+      , Math.Root.Isolation.Bisection
+      , Math.Root.Isolation.GaussSeidel
+      , Math.Root.Isolation.Narrowing
+      , Math.Root.Isolation.Core
       , Debug.Utils
 
     other-modules:

brush-strokes/src/cusps/bench/Main.hs

@@ -121,7 +121,9 @@ benchTestCase ( TestCase { testName, testBrushStroke, testCuspOptions, testStart
   before <- getMonotonicTime
   let ( _, testStrokeFnI ) = brushStrokeFunctions testBrushStroke
       ( dunno, sols ) =
-        foldMap ( \ ( i, ( _trees, ( mbCusps, defCusps ) ) ) -> ( map ( i , ) mbCusps, map ( i, ) defCusps ) ) $
+        foldMap
+          ( \ ( i, ( _trees, DoneBoxes { doneSolBoxes = defCusps, doneGiveUpBoxes = mbCusps } ) ) ->
+              ( map ( ( i , ) . snd ) mbCusps, map ( i, ) defCusps ) ) $
           IntMap.toList $
             findCuspsIn testCuspOptions testStrokeFnI $
               IntMap.fromList

brush-strokes/src/lib/Math/Algebra/Dual.hs

@@ -51,7 +51,7 @@ deriving stock instance Functor ( D k u ) => Functor ( C k u )
 -- | @D k u v@ is the space of @k@-th order germs of functions from @u@ to @v@,
 -- represented by the algebra:
 --
--- \[ \mathbb{Z}[x_1, \ldots, x_n]/(x_1, \ldots, x_n)^(k+1) \otimes_\mathbb{Z} v \]
+-- \[ \mathbb{Z}[x_1, \ldots, x_n]/(x_1, \ldots, x_n)^{k+1} \otimes_\mathbb{Z} v \]
 --
 -- when @u@ is of dimension @n@.
 type D :: Nat -> Type -> Type -> Type

brush-strokes/src/lib/Math/Bezier/Stroke.hs

@@ -597,8 +597,8 @@ outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams
       case mbCuspOptions of
         Just cuspOptions ->
           foldMap
-            ( \ ( i, ( _trees, ( potCusps, defCusps ) ) ) ->
+            ( \ ( i, ( _trees, DoneBoxes { doneSolBoxes = defCusps, doneGiveUpBoxes = potCusps } ) ) ->
-              ( map ( i , ) potCusps, map ( i , ) defCusps )
+              ( map ( ( i , ) . fst ) potCusps, map ( i , ) defCusps )
             )
             ( IntMap.toList $ findCusps cuspOptions curvesI )
         Nothing ->
@@ -1166,7 +1166,7 @@ cuspCoords eqs ( i, box )
 findCusps
   :: RootIsolationOptions N 3
   -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-  -> IntMap ( [ ( Box N, RootIsolationTree ( Box N ) ) ], ( [ Box N ], [ Box N ] ) )
+  -> IntMap ( [ ( Box N, RootIsolationTree ( Box N ) ) ], DoneBoxes N )
 findCusps opts boxStrokeData =
   findCuspsIn opts boxStrokeData $
     IntMap.fromList
@@ -1186,7 +1186,7 @@ findCuspsIn
   :: RootIsolationOptions N 3
   -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
   -> IntMap [ Box 2 ]
-  -> IntMap ( [ ( Box N, RootIsolationTree ( Box N ) ) ], ( [ Box N ], [ Box N ] ) )
+  -> IntMap ( [ ( Box N, RootIsolationTree ( Box N ) ) ], DoneBoxes N )
 findCuspsIn opts boxStrokeData initBoxes =
   IntMap.mapWithKey ( \ i -> foldMap ( isolateRootsIn opts ( eqnPiece i ) ) ) initBoxes
     where

brush-strokes/src/lib/Math/Interval.hs

@@ -6,13 +6,15 @@
 
 module Math.Interval
   ( 𝕀(𝕀), inf, sup
-  , width
+  , width, midpoint
   , scaleInterval
   , 𝕀ℝ
+  , isCanonical
   , singleton, nonDecreasing
   , inside
   , aabb
   , extendedDivide, extendedRecip
+  , intersect, bisect
   )
   where
 
@@ -20,6 +22,7 @@ module Math.Interval
 import Prelude hiding ( Num(..), Fractional(..) )
 import Data.List
   ( nub )
+import qualified Data.List.NonEmpty as NE
 
 -- acts
 import Data.Act
@@ -35,6 +38,8 @@ import Data.Group.Generics
 
 -- brush-strokes
 import Math.Algebra.Dual
+import Math.Float.Utils
+  ( prevFP )
 import Math.Interval.Internal
   ( 𝕀(𝕀), inf, sup, scaleInterval )
 import Math.Linear
@@ -51,18 +56,62 @@ import Math.Ring
 type 𝕀ℝ n = 𝕀 ( ℝ n )
 type instance D k ( 𝕀 v ) = D k v
 
+-- | An interval reduced to a single point.
 singleton :: a -> 𝕀 a
 singleton a = 𝕀 a a
 
+-- | The width of an interval.
+--
+-- NB: assumes the endpoints are neither @NaN@ nor infinity.
 width :: AbelianGroup a => 𝕀 a -> a
 width ( 𝕀 lo hi ) = hi - lo
+{-# INLINEABLE width #-}
 
 -- | Turn a non-decreasing function into a function on intervals.
 nonDecreasing :: ( a -> b ) -> 𝕀 a -> 𝕀 b
 nonDecreasing f ( 𝕀 lo hi ) = 𝕀 ( f lo ) ( f hi )
 
+-- | Does the given value lie inside the specified interval?
 inside :: Ord a => a -> 𝕀 a -> Bool
 inside x ( 𝕀 lo hi ) = x >= lo && x <= hi
+{-# INLINEABLE inside #-}
+
+-- | Is this interval canonical, i.e. it consists of either 1 or 2 floating point
+-- values only?
+isCanonical :: 𝕀 Double -> Bool
+isCanonical ( 𝕀 x_inf x_sup ) = x_inf >= prevFP x_sup
+
+-- | The midpoint of an interval.
+--
+-- NB: assumes the endpoints are neither @NaN@ nor infinity.
+midpoint :: 𝕀 Double -> Double
+midpoint ( 𝕀 x_inf x_sup ) = 0.5 * ( x_inf + x_sup )
+
+-- | Compute the intersection of two intervals.
+--
+-- Returns whether the first interval is a strict subset of the second interval
+-- (or the intersection is a single point).
+intersect :: 𝕀 Double -> 𝕀 Double -> [ ( 𝕀 Double, Bool ) ]
+intersect ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 )
+  | lo > hi
+  = [ ]
+  | otherwise
+  = [ ( 𝕀 lo hi, ( lo1 > lo2 && hi1 < hi2 ) || ( lo == hi ) ) ]
+  where
+    lo = max lo1 lo2
+    hi = min hi1 hi2
+
+-- | Bisect an interval.
+--
+-- Generally returns two sub-intervals, but returns the input interval if
+-- it is canonical (and thus bisection would be pointless).
+bisect :: 𝕀 Double -> NE.NonEmpty ( 𝕀 Double )
+bisect x@( 𝕀 x_inf x_sup )
+  | isCanonical x
+  = NE.singleton x
+  | otherwise
+  = 𝕀 x_inf x_mid NE.:| [ 𝕀 x_mid x_sup ]
+  where x_mid = midpoint x
 
 deriving via ViaAbelianGroup ( T ( 𝕀 Double ) )
   instance Semigroup ( T ( 𝕀 Double ) )

brush-strokes/src/lib/Math/Root/Isolation/Bisection.hs

@@ -0,0 +1,266 @@
+
+{-# LANGUAGE ScopedTypeVariables #-}
+
+module Math.Root.Isolation.Bisection
+  ( -- * The bisection method for root isolation
+    Bisection
+  , bisection
+
+    -- ** Configuration options
+  , BisectionCoordPicker
+  , BisectionOptions(..), defaultBisectionOptions
+
+    -- *** Helper code for picking which dimension to bisect
+  , CoordBisectionData(..)
+  , spread, normVI, maxVI
+  , sortOnArgNE
+  )
+  where
+
+-- base
+import Data.Kind
+  ( Type )
+import Data.Foldable
+ ( toList )
+import Data.Functor
+  ( (<&>) )
+import qualified Data.List.NonEmpty as NE
+  ( NonEmpty(..), cons, filter
+  , head, nonEmpty, singleton, sort
+  )
+import Data.Semigroup
+  ( Arg(..), Dual(..) )
+import GHC.TypeNats
+  ( Nat, KnownNat
+  , type (<=)
+  )
+
+-- MetaBrush
+import Math.Algebra.Dual
+  ( D )
+import Math.Interval
+import Math.Linear
+import Math.Module
+  ( Module(..) )
+import Math.Monomial
+  ( MonomialBasis(..), linearMonomial, zeroMonomial )
+import qualified Math.Ring as Ring
+import Math.Root.Isolation.Core
+
+--------------------------------------------------------------------------------
+-- Bisection
+
+-- | The bisection algorithm; see 'bisection'.
+data Bisection
+instance RootIsolationAlgorithm Bisection where
+  type instance StepDescription Bisection = ( String, Double )
+  type instance RootIsolationAlgorithmOptions Bisection n d = BisectionOptions n d
+  rootIsolationAlgorithm
+    ( BisectionOptions { canHaveSols, fallbackBisectionCoord } )
+    thisRoundHist prevRoundsHist eqs box = do
+      let ( boxes, whatBis ) =
+            bisection
+              ( canHaveSols eqs )
+              ( fallbackBisectionCoord thisRoundHist prevRoundsHist eqs )
+              box
+      return ( whatBis, boxes )
+
+-- | Options for the bisection method.
+type BisectionOptions :: Nat -> Nat -> Type
+data BisectionOptions n d =
+  BisectionOptions
+  { -- | Custom function to check whether the given box might contain solutions to the
+    -- given equations.
+    --
+    -- If you always return @True@, then we will always bisect along
+    -- the dimension picked by the 'fallbackBisectionCoord' function.
+    --
+    -- NB: only return 'False' if non-existence of solutions is guaranteed
+    -- (otherwise, the root isolation algorithm might not be consistent).
+    canHaveSols :: !( ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) -> Box n -> Bool )
+    -- | Heuristic to choose which coordinate dimension to bisect.
+    --
+    -- It's only a fallback, as we prefer to bisect along coordinate dimensions
+    -- that minimise the number of sub-boxes created.
+  , fallbackBisectionCoord :: !( BisectionCoordPicker n d )
+  }
+
+-- | A function to choose along which coordination dimension we should bisect.
+type BisectionCoordPicker n d
+  =  [ ( RootIsolationStep, Box n ) ]
+  -> BoxHistory n
+  -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
+  -> forall r. ( NE.NonEmpty ( Fin n, r ) -> ( r, String ) )
+
+-- | Default options for the bisection method.
+defaultBisectionOptions
+  :: forall n d
+  .  ( 1 <= n, BoxCt n d )
+  => Double -> Double
+  -> Box n -> BisectionOptions n d
+defaultBisectionOptions minWidth _ε_eq box =
+  BisectionOptions
+    { canHaveSols =
+        \ eqs box' ->
+          -- box(0)-consistency
+          let iRange' :: Box d
+              iRange' = eqs box' `monIndex` zeroMonomial
+          in unT ( origin @Double ) `inside` iRange'
+
+          -- box(1)-consistency
+          --let box1Options = Box1Options _ε_eq ( toList $ universe @n ) ( toList $ universe @d )
+          --in not $ null $ makeBox1Consistent _minWidth box1Options eqs box'
+
+          -- box(2)-consistency
+          --let box2Options = Box2Options _ε_eq 0.001 ( toList $ universe @n ) ( toList $ universe @d )
+          --    box'' = makeBox2Consistent _minWidth box2Options eqs box'
+          --    iRange'' :: Box d
+          --    iRange'' = eqs box'' `monIndex` zeroMonomial
+          --in unT ( origin @Double ) `inside` iRange''
+    , fallbackBisectionCoord =
+        \ _thisRoundHist _prevRoundsHist eqs possibleCoordChoices ->
+          let datPerCoord =
+                possibleCoordChoices <&> \ ( i, r ) ->
+                  CoordBisectionData
+                    { coordIndex = i
+                    , coordInterval = box `index` i
+                    , coordJacobianColumn = eqs box `monIndex` ( linearMonomial i )
+                    , coordBisectionData = r
+                    }
+
+              -- First, check if the largest dimension is over 20 times larger
+              -- than the smallest dimension; if so bisect that dimension.
+          in case sortOnArgNE ( width . coordInterval ) datPerCoord of
+                Arg _ d NE.:| [] -> ( coordBisectionData d, "" )
+                Arg w0 _ NE.:| ds ->
+                  let Arg w1 d1 = last ds
+                  in if w1 >= 20 * w0
+                     then ( coordBisectionData d1, "tooWide" )
+                     -- Otherwise, pick the coordination dimension with maximum spread
+                     -- (spread = width * Jacobian column norm).
+                     else
+                      let isTooSmall ( Arg ( Dual w ) _ ) = w < minWidth
+                      in case NE.filter ( not . isTooSmall ) $ sortOnArgNE ( Dual . spread ) datPerCoord of
+                        [] -> ( coordBisectionData d1, "tooWide'" )
+                        Arg _ d : _ -> ( coordBisectionData d, "spread" )
+              -- TODO: pick a dimension that previous Newton steps did not
+              -- manage to narrow well?
+    }
+{-# INLINEABLE defaultBisectionOptions #-}
+
+sortOnArgNE :: Ord b => ( a -> b ) -> NE.NonEmpty a -> NE.NonEmpty ( Arg b a )
+sortOnArgNE f = NE.sort . fmap ( \ a -> Arg ( f a ) a )
+{-# INLINEABLE sortOnArgNE #-}
+
+-- | A utility datatype that is useful in implementing bisection dimension picker
+-- functions ('fallbackBisectionCoord').
+type CoordBisectionData :: Nat -> Nat -> Type -> Type
+data CoordBisectionData n d r =
+  CoordBisectionData
+    { coordIndex :: !( Fin n )
+    , coordInterval :: !( 𝕀 Double )
+    , coordJacobianColumn :: !( 𝕀ℝ d )
+    , coordBisectionData :: !r
+    }
+deriving stock instance ( Show ( ℝ d ), Show r )
+                     => Show ( CoordBisectionData n d r )
+
+spread :: ( BoxCt n d, Representable Double ( ℝ d ) )
+       => CoordBisectionData n d r -> Double
+spread ( CoordBisectionData { coordInterval = cd, coordJacobianColumn = j_cd } )
+  = width cd * normVI j_cd
+{-# INLINEABLE spread #-}
+
+normVI :: ( Applicative ( Vec d ), Representable Double ( ℝ d ) ) => 𝕀ℝ d -> Double
+normVI ( 𝕀 los his ) =
+  sqrt $ sum ( nm1 <$> coordinates los <*> coordinates his )
+    where
+      nm1 :: Double -> Double -> Double
+      nm1 lo hi = max ( abs lo ) ( abs hi ) Ring.^ 2
+{-# INLINEABLE normVI #-}
+
+maxVI :: ( Applicative ( Vec d ), Representable Double ( ℝ d ) ) => 𝕀ℝ d -> Double
+maxVI ( 𝕀 los his ) =
+  maximum ( maxAbs <$> coordinates los <*> coordinates his )
+    where
+      maxAbs :: Double -> Double -> Double
+      maxAbs lo hi = max ( abs lo ) ( abs hi )
+{-# INLINEABLE maxVI #-}
+
+--------------------------------------------------------------------------------
+
+-- | Bisect the given box.
+--
+-- (The difficult part lies in determining along which coordinate
+-- dimension to bisect.)
+bisection
+  :: forall n
+  .  ( 1 <= n, KnownNat n, Representable Double ( ℝ n ) )
+  => ( Box n -> Bool )
+     -- ^ how to check whether a box contains solutions
+  -> ( forall r. NE.NonEmpty ( Fin n, r ) -> ( r, String ) )
+      -- ^ heuristic bisection coordinate picker
+  -> Box n
+      -- ^ the box to bisect
+  -> ( [ Box n ], ( String, Double ) )
+bisection canHaveSols pickFallbackBisCoord box =
+  case NE.nonEmpty solsList of
+    Nothing ->
+      -- We discarded dimensions along which bisection was useless
+      -- (because the interval was canonical in that dimension).
+      -- If there are no dimensions left, then don't do any bisection.
+      -- (TODO: we shouldn't really ever get here.)
+      ( [ box ], ( "noBis", 0 / 0 ) )
+    Just solsNE ->
+      case findFewestSols solsNE of
+        -- If there is a coordinate for which bisection results in no solutions,
+        -- or in fewer sub-boxes with solutions than any other coordinate choice,
+        -- pick that coordinate for bisection.
+        Arg nbSols ( ( i, ( mid, subBoxesWithSols ) ) NE.:| [] ) ->
+          ( subBoxesWithSols, ( "cd = " ++ show i ++ "(#subs=" ++ show nbSols ++ ")", mid ) )
+        -- Otherwise, fall back to the provided heuristic.
+        Arg _nbSols is ->
+          let ( ( mid, subBoxesWithSols ), why ) = pickFallbackBisCoord is
+          in ( subBoxesWithSols, ( why, mid ) )
+  where
+    solsList =
+      [ Arg ( fromIntegral $ length subBoxesWithSols ) ( i, ( mid, subBoxesWithSols ) )
+      | i <- toList $ universe @n
+      , let ( mid, subBoxes ) = bisectInCoord box i
+            subBoxesWithSols = NE.filter canHaveSols subBoxes
+        -- discard coordinate dimensions in which the box is a singleton
+      , length subBoxes >= 2 || null subBoxesWithSols
+      ]
+{-# INLINEABLE bisection #-}
+
+-- | Bisect a box in the given coordinate dimension.
+bisectInCoord
+  :: Representable Double ( ℝ n )
+  => Box n -> Fin n -> ( Double, NE.NonEmpty ( Box n ) )
+bisectInCoord box i =
+  let z = box `index` i
+      zs' = bisect z
+  in ( sup ( NE.head zs' )
+     , fmap ( \ z' -> set i z' box ) zs' )
+{-# INLINEABLE bisectInCoord #-}
+
+-- | Return the elements with the least argument.
+--
+-- NB: this function shortcuts as soon as it finds an element with argument 0.
+findFewestSols :: forall a. NE.NonEmpty ( Arg Word a ) -> Arg Word ( NE.NonEmpty a )
+findFewestSols ( Arg nbSols arg NE.:| args )
+  | nbSols == 0
+  = Arg 0 $ NE.singleton arg
+  | otherwise
+  = go nbSols ( NE.singleton arg ) args
+  where
+    go :: Word -> NE.NonEmpty a -> [ Arg Word a ] -> Arg Word ( NE.NonEmpty a )
+    go bestNbSolsSoFar bestSoFar [] = Arg bestNbSolsSoFar bestSoFar
+    go bestNbSolsSoFar bestSoFar ( ( Arg nbSols' arg' ) : args' )
+      | nbSols' == 0
+      = Arg 0 $ NE.singleton arg'
+      | otherwise
+      = case compare nbSols' bestNbSolsSoFar of
+          LT -> go nbSols' ( NE.singleton arg' ) args'
+          GT -> go bestNbSolsSoFar bestSoFar args'
+          EQ -> go bestNbSolsSoFar ( arg `NE.cons` bestSoFar ) args'

brush-strokes/src/lib/Math/Root/Isolation/Core.hs

@@ -0,0 +1,265 @@
+
+{-# LANGUAGE AllowAmbiguousTypes  #-}
+{-# LANGUAGE PolyKinds            #-}
+{-# LANGUAGE ScopedTypeVariables  #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+-- | Core definitions and utilities common to root isolation methods.
+module Math.Root.Isolation.Core
+  ( -- * Root isolation types
+    Box
+  , DoneBoxes(..), noDoneBoxes
+
+    -- * General typeclass for root isolation methods
+  , BoxCt
+  , RootIsolationAlgorithm(..)
+  , RootIsolationAlgorithmWithOptions(..)
+
+    -- ** Inspecting history
+  , RootIsolationStep(IsolationStep, ..)
+  , BoxHistory
+
+    -- ** Visualising history
+  , RootIsolationTree(..)
+  , showRootIsolationTree
+
+    -- * Utility functions
+  , pipeFunctionsWhileTrue
+  , forEachCoord
+  ) where
+
+-- base
+import Data.Foldable
+  ( toList )
+import Data.Kind
+  ( Type, Constraint )
+import qualified Data.List.NonEmpty as NE
+  ( NonEmpty )
+import Data.Type.Equality
+  ( (:~~:)(HRefl) )
+import Data.Typeable
+  ( Typeable, heqT )
+import GHC.TypeNats
+  ( Nat, KnownNat, type (<=) )
+import Numeric
+  ( showFFloat )
+
+-- containers
+import Data.Tree
+  ( Tree(..) )
+
+-- transformers
+import Control.Monad.Trans.State.Strict as State
+  ( State, get, put )
+import Control.Monad.Trans.Writer.CPS
+  ( Writer )
+
+-- MetaBrush
+import Math.Algebra.Dual
+  ( D )
+import Math.Interval
+import Math.Linear
+import Math.Module
+  ( Module(..) )
+import Math.Monomial
+  ( MonomialBasis(..), Deg, Vars )
+
+--------------------------------------------------------------------------------
+
+-- | An axis-aligned box in @n@-dimensions.
+type Box n = 𝕀ℝ n
+
+-- | Dimension constraints for root isolation in a system of equations:
+--
+--   - @n@: number of variables
+--   - @d@: number of equations
+--
+-- NB: we require n <= d (no support for under-constrained systems).
+--
+-- NB: in practice, this constraint should specialise away.
+type BoxCt n d =
+  ( KnownNat n, KnownNat d
+  , 1 <= n, 1 <= d, n <= d
+
+  , Show ( 𝕀ℝ n ), Show ( ℝ n )
+  , Eq ( ℝ n )
+  , Representable Double ( ℝ n )
+  , MonomialBasis ( D 1 ( ℝ n ) )
+  , Deg ( D 1 ( ℝ n ) ) ~ 1
+  , Vars ( D 1 ( ℝ n ) ) ~ n
+  , Module Double ( T ( ℝ n ) )
+  , Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) )
+  , Applicative ( Vec n )
+
+  , Ord ( ℝ d )
+  , Module Double ( T ( ℝ d ) )
+  , Representable Double ( ℝ d )
+  )
+
+-- | Boxes we are done with and will not continue processing.
+data DoneBoxes n =
+  DoneBoxes
+    { -- | Boxes which definitely contain a unique solution.
+      doneSolBoxes    :: ![ Box n ]
+      -- | Boxes which may or may not contain solutions,
+      -- and that we have stopped processing for some reason.
+    , doneGiveUpBoxes :: ![ ( Box n, String ) ]
+    }
+deriving stock instance Show ( Box n ) => Show ( DoneBoxes n )
+
+instance Semigroup ( DoneBoxes n ) where
+  DoneBoxes a1 b1 <> DoneBoxes a2 b2 = DoneBoxes ( a1 <> a2 ) ( b1 <> b2 )
+instance Monoid ( DoneBoxes n ) where
+  mempty = noDoneBoxes
+noDoneBoxes :: DoneBoxes n
+noDoneBoxes = DoneBoxes [] []
+
+--------------------------------------------------------------------------------
+-- Class for all root isolation algorithms.
+--
+-- This keeps the implementation open-ended, and allows inspection of
+-- other root isolation methods, so that heuristics can look at
+-- what happened in previous steps to decide what to do.
+
+-- | Existential wrapper over any root isolation algorithm,
+-- with the options necessary to run it.
+data RootIsolationAlgorithmWithOptions n d where
+  AlgoWithOptions
+    :: RootIsolationAlgorithm ty
+    => RootIsolationAlgorithmOptions ty n d
+    -> RootIsolationAlgorithmWithOptions n d
+
+-- | Type-class for root isolation algorithms.
+--
+-- This design keeps the set of root isolation algorithms open-ended,
+-- while retaining the ability to inspect previous steps (using the
+-- 'IsolationStep' pattern).
+type RootIsolationAlgorithm :: forall {k}. k -> Constraint
+class ( Typeable ty, Show ( StepDescription ty ) )
+    => RootIsolationAlgorithm ty where
+  -- | The type of additional information about an algorithm step.
+  --
+  -- Only really useful for debugging; gets stored in 'RootIsolationTree's.
+  type StepDescription ty
+  -- | Configuration options expected by this root isolation method.
+  type RootIsolationAlgorithmOptions ty (n :: Nat) (d :: Nat) = r | r -> ty n d
+  -- | Run one step of the root isolation method.
+  --
+  -- This gets given the equations and a box, and should attempt to
+  -- shrink the box in some way, returning smaller boxes.
+  --
+  -- Should return:
+  --
+  --  - a description of the step taken (see 'StepDescription'),
+  --  - new boxes to process (the return value of type @['Box' n]@),
+  --    which can be empty if the algorithm can prove that the input
+  --    bix does not contain any solutions;
+  --  - (as a writer side-effect) boxes to definitely stop processing; see 'DoneBoxes'.
+  rootIsolationAlgorithm
+    :: forall (n :: Nat) (d :: Nat)
+    .  BoxCt n d
+    => RootIsolationAlgorithmOptions ty n d
+        -- ^ options for this root isolation algorithm
+    -> [ ( RootIsolationStep, Box n ) ]
+        -- ^ history of the current round
+    -> BoxHistory n
+        -- ^ previous rounds history
+    -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
+      -- ^ equations
+    -> Box n
+       -- ^ box
+    -> Writer ( DoneBoxes n ) ( StepDescription ty, [ Box n ] )
+
+-- | Match on an unknown root isolation algorithm step with a known algorithm.
+pattern IsolationStep
+  :: forall (ty :: Type)
+  .  RootIsolationAlgorithm ty
+  => StepDescription ty -> RootIsolationStep
+pattern IsolationStep stepDescr
+          <- ( rootIsolationAlgorithmStep_maybe @ty -> Just stepDescr )
+  where
+    IsolationStep stepDescr = SomeRootIsolationStep @ty stepDescr
+
+-- | Helper function used to define the 'IsolationStep' pattern.
+--
+-- Inspects whether an existential 'RootIsolationStep' packs a step for
+-- the given algorithm.
+rootIsolationAlgorithmStep_maybe
+  :: forall ty. RootIsolationAlgorithm ty
+  => RootIsolationStep -> Maybe ( StepDescription ty )
+rootIsolationAlgorithmStep_maybe ( SomeRootIsolationStep @existential descr )
+  | Just HRefl <- heqT @existential @ty
+  = Just descr
+  | otherwise
+  = Nothing
+{-# INLINEABLE rootIsolationAlgorithmStep_maybe #-}
+
+-- | History for a given box: what was the outcome of previous root isolation
+-- methods?
+type BoxHistory n = [ NE.NonEmpty ( RootIsolationStep, Box n ) ]
+
+-- | A description of a step taken when isolating roots.
+data RootIsolationStep where
+  SomeRootIsolationStep
+    :: forall step
+    . ( Typeable step
+      , Show ( StepDescription step )
+      )
+    => StepDescription step
+    -> RootIsolationStep
+
+instance Show RootIsolationStep where
+  showsPrec p ( SomeRootIsolationStep stepDescr ) = showsPrec p stepDescr
+
+--------------------------------------------------------------------------------
+-- Trees recording steps taken by the algorithm, for visualisation & debugging.
+
+-- | A tree recording the steps taken when isolating roots.
+data RootIsolationTree d
+  = RootIsolationLeaf String d
+  | RootIsolationStep RootIsolationStep [ ( d, RootIsolationTree d ) ]
+
+showRootIsolationTree
+  :: ( Representable Double ( ℝ n ), Show ( Box n ) )
+  => Box n -> RootIsolationTree ( Box n ) -> Tree String
+showRootIsolationTree cand (RootIsolationLeaf why l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ why ++ " " ++ show l) []
+showRootIsolationTree cand (RootIsolationStep s ts)
+  = Node (show cand ++ " abc " ++ showArea (boxArea cand) ++ " " ++ show s) $ map (\ (c,t) -> showRootIsolationTree c t) ts
+
+boxArea :: Representable Double ( ℝ n ) => Box n -> Double
+boxArea ( 𝕀 lo hi ) =
+  product ( ( \ l h -> abs ( h - l ) ) <$> coordinates lo <*> coordinates hi )
+
+showArea :: Double -> String
+showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")"
+
+--------------------------------------------------------------------------------
+-- Utilities.
+
+-- | Run an effectful computation several times in sequence, piping its output
+-- into the next input, once for each coordinate dimension.
+forEachCoord :: forall n a m. ( KnownNat n, Monad m ) => a -> ( Fin n -> a -> m a ) -> m a
+forEachCoord a0 f = go ( toList $ universe @n ) a0
+  where
+    go [] a = return a
+    go ( i : is ) a = do
+      a' <- f i a
+      go is a'
+{-# INLINEABLE forEachCoord #-}
+
+-- | Apply each function in turn, piping the output of one function into
+-- the next.
+--
+-- Once all the functions have been applied, check whether the Bool is True.
+-- If it is, go around again with all the functions; otherwise, stop.
+pipeFunctionsWhileTrue :: [ a -> State Bool [ a ] ] -> a -> State Bool [ a ]
+pipeFunctionsWhileTrue fns = go fns
+  where
+    go [] x = do
+      doAnotherRound <- State.get
+      if doAnotherRound
+      then do { State.put False ; go fns x }
+      else return [ x ]
+    go ( f : fs ) x = do
+      xs <- f x
+      concat <$> traverse ( go fs ) xs

brush-strokes/src/lib/Math/Root/Isolation/GaussSeidel.hs

@@ -0,0 +1,239 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+
+module Math.Root.Isolation.GaussSeidel
+  ( -- * The interval Newton method with Gauss–Seidel step
+    GaussSeidel
+  , intervalGaussSeidel
+
+    -- ** Configuration options
+  , GaussSeidelOptions(..), Preconditioner(..)
+  , defaultGaussSeidelOptions
+  )
+  where
+
+-- base
+import Control.Arrow
+  ( first )
+import Data.Bifunctor
+  ( Bifunctor(bimap) )
+import Data.Coerce
+  ( coerce )
+import Data.Kind
+  ( Type )
+import Data.Foldable
+ ( toList )
+import Data.List
+  ( partition )
+import Data.Proxy
+  ( Proxy(..) )
+import Data.Type.Ord
+  ( OrderingI(..) )
+import GHC.TypeNats
+  ( Nat, KnownNat, type (<=), cmpNat )
+
+-- eigen
+import qualified Eigen.Matrix as Eigen
+  ( Matrix
+  , determinant, generate, inverse, unsafeCoeff
+  )
+
+-- transformers
+import Control.Monad.Trans.Writer.CPS
+  ( Writer, tell )
+
+-- MetaBrush
+import Math.Algebra.Dual
+  ( D )
+import Math.Epsilon
+  ( nearZero )
+import Math.Interval
+import Math.Linear
+import Math.Module
+  ( Module(..) )
+import Math.Monomial
+  ( MonomialBasis(..), linearMonomial, zeroMonomial )
+import Math.Root.Isolation.Core
+
+--------------------------------------------------------------------------------
+-- Gauss–Seidel
+
+-- | The interval Newton method with a Gauss–Seidel step; see 'intervalGaussSeidel'.
+data GaussSeidel
+instance RootIsolationAlgorithm GaussSeidel where
+  type instance StepDescription GaussSeidel = ()
+  type instance RootIsolationAlgorithmOptions GaussSeidel n d = GaussSeidelOptions n d
+  rootIsolationAlgorithm opts _thisRoundHist _prevRoundsHist eqs box = do
+    res <- intervalGaussSeidel opts eqs box
+    return ( (), res )
+  {-# INLINEABLE rootIsolationAlgorithm #-}
+
+-- | Options for the interval Gauss–Seidel method.
+type GaussSeidelOptions :: Nat -> Nat -> Type
+data GaussSeidelOptions n d =
+  GaussSeidelOptions
+    { -- | Which preconditioner to user?
+      gsPreconditioner :: !Preconditioner
+      -- | Function that projects over the equations we will consider
+      -- (the identity for a well-determined problem, or a projection for
+      -- an overdetermined system).
+    , gsPickEqs :: ( 𝕀ℝ d -> 𝕀ℝ n ) }
+
+-- | Default options for the interval Gauss–Seidel method.
+defaultGaussSeidelOptions
+  :: forall n d
+  .  ( KnownNat n, KnownNat d
+     , 1 <= n, 1 <= d, n <= d
+     , Representable Double ( ℝ n )
+     , Representable Double ( ℝ d )
+     )
+  => BoxHistory n
+  -> GaussSeidelOptions n d
+defaultGaussSeidelOptions history =
+  GaussSeidelOptions
+    { gsPreconditioner = InverseMidJacobian
+    , gsPickEqs =
+        case cmpNat @n @d Proxy Proxy of
+          EQI -> id
+          LTI ->
+            -- If there are more equations (d) than variables (n),
+            -- pick a size n subset of the variables
+            -- (go through all combinations cyclically).
+            let choices :: [ Vec n ( Fin d ) ]
+                choices = choose @d @n
+                choice :: Vec n ( Fin d )
+                choice = choices !! ( length history `mod` length choices )
+            in \ u -> tabulate \ i -> index u ( choice ! i )
+    }
+{-# INLINEABLE defaultGaussSeidelOptions #-}
+
+-- | Preconditioner to use with the interval Gauss–Seidel method.
+data Preconditioner
+  = NoPreconditioning
+  | InverseMidJacobian
+  deriving stock Show
+
+-- | Interval Newton method with Gauss–Seidel step.
+intervalGaussSeidel
+  :: forall n d
+  .  BoxCt n d
+  => GaussSeidelOptions n d
+  -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
+      -- ^ equations
+  -> 𝕀ℝ n
+     -- ^ box
+  -> Writer ( DoneBoxes n ) [ 𝕀ℝ n ]
+intervalGaussSeidel
+  ( GaussSeidelOptions { gsPreconditioner = precondMeth, gsPickEqs = pickEqs } )
+  eqs
+  box
+  | let boxMid = singleton $ boxMidpoint box
+        f' :: Vec n ( 𝕀ℝ n )
+        f' = fmap ( \ i -> pickEqs $ eqs box `monIndex` linearMonomial i ) ( universe @n )
+        f_mid = pickEqs $ eqs boxMid `monIndex` zeroMonomial
+
+  = let -- Interval Newton method: take one Gauss–Seidel step
+        -- for the system of equations f'(x) ( x - x_mid ) = - f(x_mid).
+        ( a, b ) = precondition precondMeth
+                      ( fmap boxMidpoint f' )
+                      f' ( singleton $ unT $ -1 *^ T ( boxMidpoint f_mid ) )
+
+        -- NB: we have to change coordinates, putting the midpoint of the box
+        -- at the origin, in order to take a Gauss–Seidel step.
+        gsGuesses = map ( first ( \ box' -> unT $ box' ^+^ T boxMid ) )
+                  $ gaussSeidelStep a b ( T box ^-^ T boxMid )
+    in
+         -- If the Gauss–Seidel step was a contraction, then the box
+         -- contains a unique solution (by the Banach fixed point theorem).
+         --
+         -- These boxes can thus be directly added to the solution set:
+         -- Newton's method is guaranteed to converge to the unique solution.
+         let !(done, todo) = bimap ( map fst ) ( map fst )
+                           $ partition snd gsGuesses
+         in do tell $ noDoneBoxes { doneSolBoxes = done }
+               return todo
+  where
+{-# INLINEABLE intervalGaussSeidel #-}
+
+-- | Take one interval Gauss–Seidel step for the equation \( A X = B \),
+-- refining the initial guess box for \( X \) into up to \( 2^n \) (disjoint) new boxes.
+--
+-- The boolean indicates whether the Gauss–Seidel step was a contraction.
+gaussSeidelStep
+  :: forall n
+  .  ( Representable Double ( ℝ n ), Eq ( ℝ n ) )
+  => Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \)
+  -> 𝕀ℝ n           -- ^ \( B \)
+  -> T ( 𝕀ℝ n )     -- ^ initial box \( X \)
+  -> [ ( T ( 𝕀ℝ n ), Bool ) ]
+gaussSeidelStep as b ( T x0 ) = coerce $
+  forEachCoord @n ( x0, True ) $ \ i ( x, contraction ) -> do
+    -- x_i' = ( b_i - sum { j /= i } a_ij * x_j ) / a_ii
+    x_i'0 <- ( b `index` i - sum [ ( as ! j ) `index` i * x `index` j | j <- toList ( universe @n ), j /= i ] )
+          `extendedDivide` ( ( as ! i ) `index` i )
+    ( x_i', sub_i ) <- x_i'0 `intersect` ( x `index` i )
+    return $ ( set i x_i' x, sub_i && contraction )
+-- TODO: try implementing the complete interval union Gauss–Seidel algorithm.
+-- See "Algorithm 2" in
+-- "Using interval unions to solve linear systems of equations with uncertainties"
+{-# INLINEABLE gaussSeidelStep #-}
+
+-- | The midpoint of a box.
+boxMidpoint :: Representable Double ( ℝ n ) => 𝕀ℝ n -> ℝ n
+boxMidpoint box =
+  tabulate $ \ i ->
+    let 𝕀 z_lo z_hi = index box i
+        z_mid = 0.5 * ( z_lo + z_hi )
+    in z_mid
+{-# INLINEABLE boxMidpoint #-}
+
+-- | Pre-condition the system \( AX = B \).
+precondition
+  :: forall n
+  .  ( KnownNat n, Representable Double ( ℝ n ) )
+  => Preconditioner -- ^ pre-conditioning method to use
+  -> Vec n ( ℝ n )  -- ^ entry-wise midpoint matrix of the interval Jacobian matrix
+  -> Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \)
+  -> 𝕀ℝ n           -- ^ \( B \)
+  -> ( Vec n ( 𝕀ℝ n ), 𝕀ℝ n )
+precondition meth jac_mid as b =
+  case meth of
+    NoPreconditioning
+      -> ( as, b )
+    InverseMidJacobian
+      | let mat = toEigen jac_mid
+            det = Eigen.determinant mat
+      , not $ nearZero det
+      -- (TODO: a bit wasteful to compute determinant then inverse.)
+      , let precond = Eigen.inverse mat
+            doPrecond = matMulVec ( fromEigen precond )
+      -> ( fmap doPrecond as, doPrecond b )
+      | otherwise
+      -> ( as, b )
+  where
+    toEigen :: Vec n ( ℝ n ) -> Eigen.Matrix n n Double
+    toEigen cols =
+      Eigen.generate $ \ r c ->
+        ( cols ! Fin ( fromIntegral c + 1 ) ) `index` ( Fin ( fromIntegral r + 1 ) )
+
+    fromEigen :: Eigen.Matrix n n Double -> Vec n ( ℝ n )
+    fromEigen mat =
+      fmap
+        ( \ ( Fin c ) ->
+          tabulate $ \ ( Fin r ) ->
+            Eigen.unsafeCoeff ( fromIntegral r - 1 ) ( fromIntegral c - 1 ) mat
+        )
+        ( universe @n )
+{-# INLINEABLE precondition #-}
+
+-- | Matrix multiplication \( A v \).
+matMulVec
+  :: forall n m
+  .  ( Representable Double ( ℝ n ), Representable Double ( ℝ m ) )
+  => Vec m ( ℝ n ) -- ^ columns of the matrix \( A )
+  -> 𝕀ℝ m          -- ^ vector \( v \)
+  -> 𝕀ℝ n
+matMulVec as v = tabulate $ \ r ->
+  sum [ scaleInterval ( a `index` r ) ( index v c )
+      | ( c, a ) <- toList ( (,) <$> universe @m <*> as )
+      ]
+{-# INLINEABLE matMulVec #-}

brush-strokes/src/lib/Math/Root/Isolation/Narrowing.hs

@@ -0,0 +1,490 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+
+module Math.Root.Isolation.Narrowing
+  ( -- * @box(1)@-consistency
+    Box1
+  , makeBox1Consistent
+    -- ** Configuration options
+  , Box1Options(..)
+  , defaultBox1Options
+
+  , -- *** Narrowing methods for @box(1)@-consistency
+    NarrowingMethod(..)
+  , narrowingMethods
+
+    -- **** Options for the adaptive shaving method
+  , AdaptiveShavingOptions(..)
+  , defaultAdaptiveShavingOptions
+
+    -- * @box(2)@-consistency
+  , Box2
+  , makeBox2Consistent
+
+    -- ** Configuration options
+  , Box2Options(..)
+  , defaultBox2Options
+  )
+  where
+
+-- base
+import Control.Monad
+  ( when )
+import Data.Foldable
+  ( toList )
+import qualified Data.List.NonEmpty as NE
+  ( toList )
+import GHC.TypeNats
+  ( KnownNat )
+
+-- transformers
+import Control.Monad.Trans.State.Strict as State
+  ( State, evalState, get, put )
+
+-- brush-strokes
+import Math.Algebra.Dual
+  ( D )
+import Math.Float.Utils
+  ( succFP, prevFP )
+import Math.Interval
+import Math.Linear
+  ( ℝ
+  , Representable
+  , Fin, set, index, universe
+  )
+import Math.Monomial
+  ( MonomialBasis(..), Deg, Vars
+  , zeroMonomial, linearMonomial
+  )
+import Math.Root.Isolation.Core
+
+--------------------------------------------------------------------------------
+-- Box-consistency driver code
+
+-- | A @box(1)@-consistency enforcing algorithm; see 'makeBox1Consistent'.
+data Box1
+instance RootIsolationAlgorithm Box1 where
+  type instance StepDescription Box1 = ()
+  type instance RootIsolationAlgorithmOptions Box1 n d = Box1Options n d
+  rootIsolationAlgorithm opts _thisRoundHist _prevRoundsHist eqs box =
+    return $ ( (), makeBox1Consistent opts eqs box )
+  {-# INLINEABLE rootIsolationAlgorithm #-}
+
+-- | A @box(2)@-consistency enforcing algorithm; see 'makeBox1Consistent'.
+data Box2
+instance RootIsolationAlgorithm Box2 where
+  type instance StepDescription Box2 = ()
+  type instance RootIsolationAlgorithmOptions Box2 n d = Box2Options n d
+  rootIsolationAlgorithm opts _thisRoundHist _prevRoundsHist eqs box =
+    return ( () , [ makeBox2Consistent opts eqs box ] )
+  {-# INLINEABLE rootIsolationAlgorithm #-}
+
+-- | Options for the @box(1)@-consistency method.
+data Box1Options n d =
+  Box1Options
+    { box1EpsEq           :: !Double
+    , box1EpsBis          :: !Double
+    , box1CoordsToNarrow  :: !( [ Fin n ] )
+    , box1EqsToUse        :: !( [ Fin d ] )
+    , box1NarrowingMethod :: !NarrowingMethod
+    }
+  deriving stock Show
+
+-- | Options for the @box(2)@-consistency method.
+data Box2Options n d =
+  Box2Options
+    { box2Box1Options :: !( Box1Options n d )
+    , box2EpsEq       :: !Double
+    , box2LambdaMin   :: !Double
+    }
+  deriving stock Show
+
+-- | Default options for the @box(1)@-consistency method.
+defaultBox1Options
+  :: forall n d
+  .  ( KnownNat n, KnownNat d )
+  => Double -- ^ minimum width of boxes (don't bisect further)
+  -> Double -- ^ threshold for progress
+  -> Box1Options n d
+defaultBox1Options minWidth ε_eq =
+  Box1Options
+    { box1EpsEq           = ε_eq
+    , box1EpsBis          = minWidth
+    , box1CoordsToNarrow  = toList $ universe @n
+    , box1EqsToUse        = toList $ universe @d
+    , box1NarrowingMethod = Kubica
+    }
+{-# INLINEABLE defaultBox1Options #-}
+
+-- | Default options for the @box(2)@-consistency method.
+defaultBox2Options
+  :: forall n d
+  .  ( KnownNat n, KnownNat d )
+  => Double -- ^ minimum width of boxes (don't bisect further)
+  -> Double -- ^ threshold for progress
+  -> Box2Options n d
+defaultBox2Options minWidth ε_eq =
+  Box2Options
+    { box2Box1Options = defaultBox1Options minWidth ε_eq
+    , box2EpsEq       = ε_eq
+    , box2LambdaMin   = 0.001
+    }
+
+-- | An implementation of "bc_enforce" from the paper
+-- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
+--
+-- See also
+-- "Presentation of a highly tuned multithreaded interval solver for underdetermined and well-determined nonlinear systems"
+makeBox1Consistent
+  :: ( KnownNat n, Representable Double ( ℝ n )
+     , MonomialBasis ( D 1 ( ℝ n ) )
+     , Deg ( D 1 ( ℝ n ) ) ~ 1
+     , Vars ( D 1 ( ℝ n ) ) ~ n
+     , Representable Double ( ℝ d )
+     )
+  => Box1Options n d
+  -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
+  -> 𝕀ℝ n -> [ 𝕀ℝ n ]
+makeBox1Consistent box1Options eqs x =
+  ( `State.evalState` False ) $
+    pipeFunctionsWhileTrue ( allNarrowingOperators box1Options eqs ) x
+{-# INLINEABLE defaultBox2Options #-}
+
+-- | An implementation of "bound-consistency" from the paper
+-- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
+makeBox2Consistent
+  :: forall n d
+  .  ( KnownNat n, Representable Double ( ℝ n )
+     , MonomialBasis ( D 1 ( ℝ n ) )
+     , Deg ( D 1 ( ℝ n ) ) ~ 1
+     , Vars ( D 1 ( ℝ n ) ) ~ n
+     , Representable Double ( ℝ d )
+     )
+  => Box2Options n d
+  -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
+  -> 𝕀ℝ n -> 𝕀ℝ n
+makeBox2Consistent (Box2Options box1Options ε_eq λMin) eqs x0
+  = ( `State.evalState` False ) $ doLoop 0.25 x0
+  where
+    doBox1 :: 𝕀ℝ n -> [ 𝕀ℝ n ]
+    doBox1 = makeBox1Consistent box1Options eqs
+    doLoop :: Double -> 𝕀ℝ n -> State Bool ( 𝕀ℝ n )
+    doLoop λ x = do
+      x'' <- forEachCoord @n x $ boundConsistency λ
+      modified <- State.get
+      let λ' = if modified then λ else 0.5 * λ
+      if λ' < λMin
+      then return x''
+      else do { State.put False ; doLoop λ' x'' }
+
+    boundConsistency :: Double -> Fin n -> 𝕀ℝ n -> State Bool ( 𝕀ℝ n )
+    boundConsistency λ i box = do
+      let x@( 𝕀 x_inf x_sup ) = getter box
+          c1 = ( 1 - λ ) * x_inf + λ * x_sup
+          c2 = λ * x_inf + ( 1 - λ ) * x_sup
+          x'_inf =
+            case doBox1 ( setter ( 𝕀 x_inf c1 ) box ) of
+              []  -> c1
+              x's -> minimum $ map ( inf . getter ) x's
+          x'_sup =
+            case doBox1 ( setter ( 𝕀 c2 x_sup ) box ) of
+              []  -> c2
+              x's -> maximum $ map ( sup . getter ) x's
+          x' = 𝕀 x'_inf x'_sup
+      when ( width x - width x' >= ε_eq ) $
+        State.put True
+      return $ setter x' box
+        where
+          getter = ( `index` i )
+          setter = set i
+
+--------------------------------------------------------------------------------
+-- Narrowing methods
+
+-- | The narrowing method to use to enforce @box(1)@-consistency.
+data NarrowingMethod
+  -- | Algorithm 5 from the paper
+  -- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
+  -- (Bartłomiej Jacek Kubica, 2017)
+  = Kubica
+  -- | Two-sided shaving @sbc@ from the paper
+  -- "A Data-Parallel Algorithm to Reliably Solve Systems of Nonlinear Equations",
+  -- (Goldsztejn & Goualard, 2008)
+  | TwoSidedShaving
+  -- | @sbc3ag@ (adaptive guessing) shaving, from the paper
+  -- "Box Consistency through Adaptive Shaving"
+  -- (Goldsztejn & Goualard, 2010).
+  | AdaptiveShaving AdaptiveShavingOptions
+  deriving stock Show
+
+narrowingMethods
+  :: Double -> Double
+  -> NarrowingMethod
+  -> [ ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) -> 𝕀 Double -> [ 𝕀 Double ] ]
+narrowingMethods ε_eq ε_bis Kubica
+  = [ leftNarrow ε_eq ε_bis, rightNarrow ε_eq ε_bis ]
+narrowingMethods ε_eq ε_bis ( AdaptiveShaving opts )
+  = [ leftShave ε_eq ε_bis opts, rightNarrow ε_eq ε_bis ]
+    -- TODO: haven't implemented right shaving yet
+narrowingMethods _ε_eq ε_bis TwoSidedShaving
+  = [ sbc ε_bis ]
+{-# INLINE narrowingMethods #-}
+
+allNarrowingOperators
+  :: forall n d
+  .  ( KnownNat n, Representable Double ( ℝ n )
+     , MonomialBasis ( D 1 ( ℝ n ) )
+     , Deg ( D 1 ( ℝ n ) ) ~ 1
+     , Vars ( D 1 ( ℝ n ) ) ~ n
+     , Representable Double ( ℝ d )
+     )
+  => Box1Options n d
+  -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
+  -> [ 𝕀ℝ n -> State Bool [ 𝕀ℝ n ] ]
+allNarrowingOperators ( Box1Options ε_eq ε_bis coordsToNarrow eqsToUse narrowingMethod ) eqs =
+  [ \ cand ->
+    let getter = ( `index` coordIndex )
+        setter = set coordIndex
+        newCands = map ( `setter` cand )
+                 $ narrowFn ( \ box -> ff' coordIndex eqnIndex $ setter box cand )
+                     ( getter cand )
+    in do
+      -- Record when we achieved a meaningful reduction,
+      -- so that we continue trying further narrowings.
+      when ( ( width ( getter cand ) - sum ( map ( width . getter ) newCands ) ) >= ε_eq ) $
+        -- NB: making this return 'True' less often seems slightly beneficial?
+        -- Further investigation needed.
+        State.put True
+      return newCands
+  | narrowFn <- narrowingMethods ε_eq ε_bis narrowingMethod
+  , coordIndex <- coordsToNarrow
+  , eqnIndex <- eqsToUse
+  ]
+  where
+    ff' :: Fin n -> Fin d -> 𝕀ℝ n -> ( 𝕀 Double, 𝕀 Double )
+    ff' i d ts =
+      let df = eqs ts
+          f, f' :: 𝕀ℝ d
+          f = df `monIndex` zeroMonomial
+          f' = df `monIndex` linearMonomial i
+      in ( f `index` d, f' `index` d )
+{-# INLINEABLE allNarrowingOperators #-}
+
+--------------------------------------------------------------------------------
+-- Kubica's algorithm.
+
+-- Use the univariate interval Newton method to narrow from the left
+-- a candidate interval.
+--
+-- See Algorithm 5 (Procedure left_narrow) in
+-- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
+-- by Bartłomiej Jacek Kubica, 2017
+leftNarrow :: Double
+           -> Double
+           -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) )
+           -> 𝕀 Double
+           -> [ 𝕀 Double ]
+leftNarrow ε_eq ε_bis ff' = left_narrow
+  where
+    left_narrow ( 𝕀 x_inf x_sup ) =
+      go x_sup ( 𝕀 x_inf ( if x_inf == x_sup then x_inf else succFP x_inf ) )
+    go x_sup x_left =
+      let ( f_x_left, _f'_x_left ) = ff' x_left
+      in
+        if inf f_x_left <= 0 && sup f_x_left >= 0
+        then [ 𝕀 ( inf x_left ) x_sup ]
+        else
+          let x = 𝕀 ( sup x_left ) x_sup
+              ( _f_x, f'_x ) = ff' x
+              x's = do δ <- f_x_left `extendedDivide` f'_x
+                       let x_new = x_left - δ
+                       map fst $ x_new `intersect` x
+          in
+            if | null x's
+               -> []
+               | ( width x - sum ( map width x's ) ) < ε_eq
+               -> x's
+               | otherwise
+               -> do
+                x' <- x's
+                if sup x' - inf x' < ε_bis
+                then return x'
+                else left_narrow =<< NE.toList ( bisect x' )
+
+-- TODO: de-duplicate with 'leftNarrow'?
+rightNarrow :: Double
+            -> Double
+            -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) )
+            -> 𝕀 Double
+            -> [ 𝕀 Double ]
+rightNarrow ε_eq ε_bis ff' = right_narrow
+  where
+    right_narrow ( 𝕀 x_inf x_sup ) =
+      go x_inf ( 𝕀 ( if x_inf == x_sup then x_sup else prevFP x_sup ) x_sup )
+    go x_inf x_right =
+      let ( f_x_right, _f'_x_right ) = ff' x_right
+      in
+        if inf f_x_right <= 0 && sup f_x_right >= 0
+        then [ 𝕀 x_inf ( sup x_right ) ]
+        else
+          let x = 𝕀 x_inf ( inf x_right )
+              ( _f_x, f'_x ) = ff' x
+              x's = do δ <- f_x_right `extendedDivide` f'_x
+                       let x_new = x_right - δ
+                       map fst $ x_new `intersect` x
+          in
+            if | null x's
+               -> []
+               | ( width x - sum ( map width x's ) ) < ε_eq
+               -> x's
+               | otherwise
+               -> do
+                x' <- x's
+                if sup x' - inf x' < ε_bis
+                then return x'
+                else right_narrow =<< NE.toList ( bisect x' )
+
+--------------------------------------------------------------------------------
+-- Adaptive shaving.
+
+data AdaptiveShavingOptions
+  = AdaptiveShavingOptions
+  { γ_init, σ_good, σ_bad, β_good, β_bad :: !Double
+  }
+  deriving stock Show
+
+defaultAdaptiveShavingOptions :: AdaptiveShavingOptions
+defaultAdaptiveShavingOptions =
+  AdaptiveShavingOptions
+    { γ_init = 0.25
+    , σ_good = 0.25
+    , σ_bad  = 0.75
+    , β_good = 1.5
+    , β_bad  = 0.7
+    }
+
+-- | Algorithm @lnar_sbc3ag@ (adaptive guessing) from the paper
+--    "Box Consistency through Adaptive Shaving" (Goldsztejn & Goualard, 2010).
+leftShave :: Double
+          -> Double
+          -> AdaptiveShavingOptions
+          -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) )
+          -> 𝕀 Double
+          -> [ 𝕀 Double ]
+leftShave ε_eq ε_bis
+  ( AdaptiveShavingOptions { γ_init, σ_good, σ_bad, β_good, β_bad } )
+  ff' i0 =
+    left_narrow γ_init i0
+  where
+    w0 = width i0
+    left_narrow :: Double -> 𝕀 Double -> [ 𝕀 Double ]
+    left_narrow γ i@( 𝕀 x_inf x_sup )
+      -- Stop if the box is too small.
+      | width i < ε_bis
+      = [ i ]
+      | otherwise
+      = go γ x_sup ( 𝕀 x_inf ( if x_inf == x_sup then x_inf else succFP x_inf ) )
+    go :: Double -> Double -> 𝕀 Double -> [ 𝕀 Double ]
+    go γ x_sup x_left =
+      let ( f_x_left, _f'_x_left ) = ff' x_left
+      in
+        if 0 `inside` f_x_left
+        -- Box-consistency achieved; finish.
+        then [ 𝕀 ( inf x_left ) x_sup ]
+        else
+          -- Otherwise, try to shave off a chunk on the left of the interval.
+          let x' = 𝕀 ( sup x_left ) x_sup
+              inf_guess = sup x_left
+              sup_guess = min x_sup ( inf_guess + γ * w0 ) -- * width x' )
+              guess = 𝕀 inf_guess sup_guess
+                -- NB: this always uses the initial width (to "avoid asymptotic behaviour" according to the paper)
+              ( f_guess, f'_guess ) = ff' guess
+              x_minus_guess = 𝕀 ( min x_sup ( succFP $ sup guess ) ) x_sup
+          in if not ( 0 `inside` f_guess )
+             then
+               -- We successfully shaved "guess" off; go round again after removing it.
+               -- TODO: here we could go back to the top with a new "w" maybe?
+               left_narrow ( β_good * γ ) x_minus_guess
+             else
+                -- Do a Newton step to try to reduce the guess interval.
+                -- Starting from "guess", we get a collection of sub-intervals
+                -- "guesses'" refining where the function can be zero.
+                let guesses' :: [ ( 𝕀 Double ) ]
+                    guesses' = do
+                      δ <- f_x_left `extendedDivide` f'_guess
+                      let guess' = singleton ( inf guess ) - δ
+                      map fst $ guess' `intersect` guess
+                    w_guess = width guess
+                    w_guesses'
+                      | null guesses'
+                      = 0
+                      | otherwise
+                      = sup_guess - minimum ( map inf guesses' )
+                    γ' | w_guesses' < σ_good * w_guess
+                       -- Good improvement, try larger guesses in the future.
+                       = β_good * γ
+                       | w_guesses' > σ_bad * w_guess
+                       -- Poor improvement, try smaller guesses in the future.
+                       = β_bad * γ
+                       | otherwise
+                       -- Otherwise, keep the γ factor the same.
+                       = γ
+                    xs' = x_minus_guess : guesses'
+                in if ( width x' - sum ( map width xs' ) ) < ε_eq
+                   then xs'
+                   else left_narrow γ' =<< xs'
+{-# INLINEABLE leftShave #-}
+
+--------------------------------------------------------------------------------
+-- Two-sided shaving.
+
+-- | @sbc@ algorithm from the paper
+--
+-- "A Data-Parallel Algorithm to Reliably Solve Systems of Nonlinear Equations",
+-- (Frédéric Goualard, Alexandre Goldsztejn, 2008).
+sbc :: Double
+    -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) )
+    -> 𝕀 Double -> [ 𝕀 Double ]
+sbc ε ff' = go
+  where
+    go :: 𝕀 Double -> [ 𝕀 Double ]
+    go x@( 𝕀 x_l x_r )
+      | width x <= ε
+      = [ x ]
+      | otherwise
+      = let x_mid = 0.5 * ( x_l + x_r )
+            ( left_done, left_todo )
+              | 0 `inside` ( fst $ ff' ( 𝕀 x_l x_lp ) )
+              = ( True, [ ] )
+              | not $ 0 `inside` ( fst $ ff' i_l )
+              = ( False, [ ] )
+              | otherwise
+              = let
+                   xls = do
+                      let l = 𝕀 x_lp x_lp
+                      δ <- fst ( ff' l ) `extendedDivide` snd ( ff' i_l )
+                      map fst $ ( l - δ ) `intersect` i_l
+                in ( False, xls )
+              where x_lp = min ( succFP x_l ) x_mid
+                    i_l = 𝕀 x_lp x_mid
+            ( right_done, right_todo )
+              | 0 `inside` ( fst $ ff' ( 𝕀 x_rm x_r ) )
+              = ( True, [ ] )
+              | not $ 0 `inside` ( fst $ ff' i_r )
+              = ( False, [ ] )
+              | otherwise
+              = let
+                   xrs = do
+                      let r = 𝕀 x_rm x_rm
+                      δ <- fst ( ff' r ) `extendedDivide` snd ( ff' i_r )
+                      map fst $ ( r - δ ) `intersect` i_r
+                in ( False, xrs )
+              where x_rm = max ( prevFP x_r ) x_mid
+                    i_r = 𝕀 x_mid x_rm
+        in do let lefts' = if left_done
+                           then [ 𝕀 x_l x_mid ]
+                           else go =<< left_todo
+                  rights' = if right_done
+                            then [ 𝕀 x_mid x_r ]
+                            else go =<< right_todo
+              lefts' ++ rights'
+{-# INLINEABLE sbc #-}