WIP: add Walter's LP approach to interval Newton

brush-strokes/brush-strokes.cabal

@@ -9,6 +9,9 @@ build-type:     Simple
 description:
   Computing brush strokes using Bézier curves.
 
+extra-source-files:
+  cbits/**/*.{cpp, hpp}
+
 flag use-fma
   description: Use fused-muliply add instructions to implement interval arithmetic.
   default: True
@@ -143,8 +146,11 @@ library
       , Math.Root.Isolation.Bisection
       , Math.Root.Isolation.Core
       , Math.Root.Isolation.Degree
-      , Math.Root.Isolation.GaussSeidel
       , Math.Root.Isolation.Narrowing
+      , Math.Root.Isolation.Newton
+      , Math.Root.Isolation.Newton.GaussSeidel
+      , Math.Root.Isolation.Newton.LP
+      , Math.Root.Isolation.Utils
       , Debug.Utils
 
     other-modules:
@@ -186,6 +192,22 @@ library
       , transformers
            >= 0.5.6.2 && < 0.7
 
+    -- Extra C++ code for 2D linear systems of interval equations
+    include-dirs:
+      cbits
+    cxx-sources:
+      cbits/lp_2d.cpp
+    cxx-options:
+      -std=c++20
+      -mavx2 -mfma
+      -frounding-math -fno-math-errno -fno-signed-zeros
+      -fno-trapping-math -fno-signaling-nans
+      -Wno-unused-result -Wsign-compare -Wno-switch
+      -march=native -mtune=native
+      -O3 -DNDEBUG
+    build-depends:
+      system-cxx-std-lib
+
 --executable convert-metafont
 --
 --    import:

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

@@ -28,10 +28,12 @@ module Math.Root.Isolation
   , BisectionOptions(..), BisectionCoordPicker
   , defaultBisectionOptions
 
-    -- * The interval Newton method with Gauss–Seidel step
-  , GaussSeidel
+    -- * The interval Newton method
+  , NewtonOptions(..)
+  , defaultNewtonOptions
+    -- ** Options for the Gauss–Seidel step
   , GaussSeidelOptions(..), Preconditioner(..)
-  , defaultGaussSeidelOptions
+  , GaussSeidelUpdateMethod(..)
 
     -- * Box-consistency methods
 
@@ -79,8 +81,9 @@ import Math.Monomial
 
 import Math.Root.Isolation.Bisection
 import Math.Root.Isolation.Core
-import Math.Root.Isolation.GaussSeidel
 import Math.Root.Isolation.Narrowing
+import Math.Root.Isolation.Newton
+import Math.Root.Isolation.Newton.GaussSeidel
 
 --------------------------------------------------------------------------------
 
@@ -134,16 +137,16 @@ defaultRootIsolationAlgorithms minWidth ε_eq history box =
     -- Currently: we try an interval Gauss–Seidel.
     -- (box(1)- and box(2)-consistency don't seem to help when using
     -- the complete interval union Gauss–Seidel step)
-    _ -> Right $ AlgoWithOptions @GaussSeidel _gsOptions
+    _ -> Right $ AlgoWithOptions @Newton _newtonOptions
             NE.:| []
 
   where
     verySmall = and $ ( \ cd -> width cd <= minWidth ) <$> coordinates box
 
-    _bisOptions = defaultBisectionOptions @n @d minWidth ε_eq box
-    _gsOptions = defaultGaussSeidelOptions @n @d history
-    _box1Options = defaultBox1Options @n @d ( minWidth * 100 ) ε_eq
-    _box2Options = ( defaultBox2Options @n @d minWidth ε_eq ) { box2LambdaMin = 0.001 }
+    _bisOptions    = defaultBisectionOptions @n @d minWidth ε_eq box
+    _newtonOptions = NewtonLP -- defaultNewtonOptions @n @d history
+    _box1Options   = defaultBox1Options @n @d ( minWidth * 100 ) ε_eq
+    _box2Options   = ( defaultBox2Options @n @d minWidth ε_eq ) { box2LambdaMin = 0.001 }
 
     -- Did we reduce the box width by at least ε_eq
     -- in at least one of the coordinates?

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

@@ -78,7 +78,8 @@ type Box n = 𝕀ℝ n
 -- NB: we require n <= d (no support for under-constrained systems).
 --
 -- NB: in practice, this constraint should specialise away.
-type BoxCt n d =
+type BoxCt n d = ( n ~ 2, d ~ 3 )
+{-
   ( KnownNat n, KnownNat d
   , 1 <= n, 1 <= d, n <= d
 
@@ -95,7 +96,7 @@ type BoxCt n d =
   , Module Double ( T ( ℝ d ) )
   , Representable Double ( ℝ d )
   )
-
+-}
 -- | Boxes we are done with and will not continue processing.
 data DoneBoxes n =
   DoneBoxes

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

@@ -0,0 +1,171 @@
+{-# LANGUAGE ScopedTypeVariables  #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+module Math.Root.Isolation.Newton
+  ( -- * The interval Newton method,
+    -- with Gauss–Seidel step or explicit linear programming
+    Newton
+  , intervalNewton
+
+    -- ** Configuration options
+  , NewtonOptions(..)
+  , defaultNewtonOptions
+  )
+  where
+
+-- base
+import Prelude hiding ( unzip )
+import Control.Arrow
+  ( first )
+import Data.Bifunctor
+  ( Bifunctor(bimap) )
+import Data.Kind
+  ( Type )
+import Data.List
+  ( partition )
+import GHC.TypeNats
+  ( Nat, KnownNat, type (<=) )
+
+-- transformers
+import Control.Monad.Trans.Writer.CPS
+  ( Writer, tell )
+
+-- MetaBrush
+import Math.Algebra.Dual
+  ( D )
+import Math.Interval
+import Math.Linear
+import Math.Module
+  ( Module(..) )
+import Math.Monomial
+  ( MonomialBasis(..), linearMonomial, zeroMonomial )
+import Math.Root.Isolation.Core
+import Math.Root.Isolation.Newton.GaussSeidel
+import Math.Root.Isolation.Newton.LP
+import Math.Root.Isolation.Utils
+
+--------------------------------------------------------------------------------
+-- Interval Newton
+
+-- | The interval Newton method; see 'intervalNewton'.
+data Newton
+instance BoxCt n d => RootIsolationAlgorithm Newton n d where
+  type instance StepDescription Newton = ()
+  type instance RootIsolationAlgorithmOptions Newton n d = NewtonOptions n d
+  rootIsolationAlgorithm opts _thisRoundHist _prevRoundsHist eqs box = do
+    res <- intervalNewton @n @d opts eqs box
+    return ( (), res )
+  {-# INLINEABLE rootIsolationAlgorithm #-}
+  {-# SPECIALISE rootIsolationAlgorithm
+        :: RootIsolationAlgorithmOptions Newton 2 3
+        -> [ ( RootIsolationStep, Box 2 ) ]
+        -> BoxHistory 2
+        -> ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) )
+        -> Box 2
+        -> Writer ( DoneBoxes 2 ) ( StepDescription Newton, [ Box 2 ] ) #-}
+    -- NB: including this to be safe. The specialiser seems to sometimes
+    -- be able to generate this specialisation on its own, and sometimes not.
+
+-- | Options for the interval Newton method.
+type NewtonOptions :: Nat -> Nat -> Type
+data NewtonOptions n d where
+  -- | Use the Gauss–Seidel method to solve linear systems.
+  NewtonGaussSeidel
+    :: GaussSeidelOptions n d -> NewtonOptions n d
+  -- | Use linear programming to solve linear systems (2 dimensions only).
+  NewtonLP
+    :: NewtonOptions 2 d
+
+-- | Default options for the interval Newton method.
+defaultNewtonOptions
+  :: forall n d
+  .  ( KnownNat n, KnownNat d
+     , 1 <= n, 1 <= d, n <= d
+     , Representable Double ( ℝ n )
+     , Representable Double ( ℝ d )
+     )
+  => BoxHistory n
+  -> NewtonOptions n d
+defaultNewtonOptions history =
+  NewtonGaussSeidel $ defaultGaussSeidelOptions history
+{-# INLINEABLE defaultNewtonOptions #-}
+
+-- | Interval Newton method with Gauss–Seidel step.
+intervalNewton
+  :: forall n d
+  .  BoxCt n d
+  => NewtonOptions n d
+  -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
+      -- ^ equations
+  -> 𝕀ℝ n
+     -- ^ box
+  -> Writer ( DoneBoxes n ) [ 𝕀ℝ n ]
+intervalNewton opts eqs x = case opts of
+  NewtonGaussSeidel
+    ( GaussSeidelOptions
+      { gsPreconditioner = precondMeth
+      , gsPickEqs = pickEqs
+      , gsUpdate
+      } ) ->
+    let x_mid = singleton $ boxMidpoint x
+        f :: 𝕀ℝ n -> 𝕀ℝ n
+        f = \ x_0 -> pickEqs $ eqs x_0 `monIndex` zeroMonomial
+        f'_x :: Vec n ( 𝕀ℝ n )
+        f'_x = fmap ( \ i -> pickEqs $ eqs x `monIndex` linearMonomial i ) ( universe @n )
+
+        -- Interval Newton method: take one Gauss–Seidel step
+        -- for the system of equations f'(x) ( x - x_mid ) = - f(x_mid).
+        minus_f_x_mid = unT $ -1 *^ T ( boxMidpoint $ f x_mid )
+
+        -- Precondition the above linear system into A ( x - x_mid ) = B.
+        ( a, b ) = precondition precondMeth
+                      f'_x ( singleton minus_f_x_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 ( \ x' -> unT $ x' ^+^ T x_mid ) )
+                  $ gaussSeidelUpdate gsUpdate a b ( T x ^-^ T x_mid )
+        ( done, todo ) = bimap ( map fst ) ( map fst )
+                       $ partition snd gsGuesses
+    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.
+      do tell $ noDoneBoxes { doneSolBoxes = done }
+         return todo
+  NewtonLP ->
+    -- TODO: reduce duplication with the above.
+    let x_mid = singleton $ boxMidpoint x
+        f :: 𝕀ℝ 2 -> 𝕀ℝ d
+        f = \ x_0 -> eqs x_0 `monIndex` zeroMonomial
+        f'_x :: Vec 2 ( 𝕀ℝ d )
+        f'_x = fmap ( \ i -> eqs x `monIndex` linearMonomial i ) ( universe @2 )
+
+        minus_f_x_mid = unT $ -1 *^ T ( boxMidpoint $ f x_mid )
+        ( a, b ) = ( f'_x, singleton minus_f_x_mid )
+        lpGuesses = map ( first ( \ x' -> unT $ x' ^+^ T x_mid ) )
+                  $ solveIntervalLinearEquations a b ( T x ^-^ T x_mid )
+        ( done, todo ) = bimap ( map fst ) ( map fst )
+                       $ partition snd lpGuesses
+    in do tell $ noDoneBoxes { doneSolBoxes = done }
+          return todo
+{-# INLINEABLE intervalNewton #-}
+{-
+
+mbDeg = topologicalDegree 0.005 f x
+det = case f'_x of
+        Vec [ c1, c2 ] ->
+          let a_11 = c1 `index` Fin 1
+              a_12 = c2 `index` Fin 1
+              a_21 = c1 `index` Fin 2
+              a_22 = c2 `index` Fin 2
+          in a_11 * a_22 - a_12 * a_21
+        _ -> error "TODO: just testing n=2 here"
+
+if | not $ 0 ∈ det
+   , mbDeg == Just 0
+   -> return []
+     -- If the Jacobian is invertible over the box, then the topological
+     -- degree tells us exactly how many solutions there are in the box.
+-}

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

@@ -1,33 +1,26 @@
-{-# LANGUAGE ScopedTypeVariables  #-}
-{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE ScopedTypeVariables #-}
 
-module Math.Root.Isolation.GaussSeidel
-  ( -- * The interval Newton method with Gauss–Seidel step
-    GaussSeidel
-  , intervalGaussSeidel
-
-    -- ** Configuration options
-  , GaussSeidelOptions(..), Preconditioner(..)
+-- | The Gauss–Seidel method for solving systems of
+-- interval linear equations
+module Math.Root.Isolation.Newton.GaussSeidel
+  ( -- *  Gauss–Seidel step
+    gaussSeidelUpdate
+    -- ** Options for the Gauss–Seidel method
+  , GaussSeidelOptions(..)
+  , Preconditioner(..), GaussSeidelUpdateMethod(..)
   , defaultGaussSeidelOptions
+  , precondition
   )
   where
 
 -- base
 import Prelude hiding ( unzip )
-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.List.NonEmpty
-  ( unzip )
+  ( toList )
 import Data.Proxy
   ( Proxy(..) )
 import Data.Type.Ord
@@ -41,44 +34,16 @@ import qualified Eigen.Matrix as Eigen
   , generate, inverse, unsafeCoeff
   )
 
--- transformers
-import Control.Monad.Trans.Writer.CPS
-  ( Writer, tell )
-
 -- MetaBrush
-import Math.Algebra.Dual
-  ( D )
 import Math.Interval
 import Math.Linear
-import Math.Module
-  ( Module(..) )
-import Math.Monomial
-  ( MonomialBasis(..), linearMonomial, zeroMonomial )
 import Math.Root.Isolation.Core
+import Math.Root.Isolation.Utils
 
 --------------------------------------------------------------------------------
 -- Gauss–Seidel
 
--- | The interval Newton method with a Gauss–Seidel step; see 'intervalGaussSeidel'.
-data GaussSeidel
-instance BoxCt n d => RootIsolationAlgorithm GaussSeidel n d where
-  type instance StepDescription GaussSeidel = ()
-  type instance RootIsolationAlgorithmOptions GaussSeidel n d = GaussSeidelOptions n d
-  rootIsolationAlgorithm opts _thisRoundHist _prevRoundsHist eqs box = do
-    res <- intervalGaussSeidel @n @d opts eqs box
-    return ( (), res )
-  {-# INLINEABLE rootIsolationAlgorithm #-}
-  {-# SPECIALISE rootIsolationAlgorithm
-        :: RootIsolationAlgorithmOptions GaussSeidel 2 3
-        -> [ ( RootIsolationStep, Box 2 ) ]
-        -> BoxHistory 2
-        -> ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) )
-        -> Box 2
-        -> Writer ( DoneBoxes 2 ) ( StepDescription GaussSeidel, [ Box 2 ] ) #-}
-    -- NB: including this to be safe. The specialiser seems to sometimes
-    -- be able to generate this specialisation on its own, and sometimes not.
-
--- | Options for the interval Gauss–Seidel method.
+-- | Options for the Gauss–Seidel method.
 type GaussSeidelOptions :: Nat -> Nat -> Type
 data GaussSeidelOptions n d =
   GaussSeidelOptions
@@ -98,7 +63,7 @@ data GaussSeidelUpdateMethod
   | GS_Complete
   deriving stock Show
 
--- | Default options for the interval Gauss–Seidel method.
+-- | Default options for the Gauss–Seidel step.
 defaultGaussSeidelOptions
   :: forall n d
   .  ( KnownNat n, KnownNat d
@@ -133,72 +98,6 @@ data Preconditioner
   | InverseMidpoint
   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
-    , gsUpdate
-    } )
-  eqs
-  x
-  | let x_mid = singleton $ boxMidpoint x
-        f :: 𝕀ℝ n -> 𝕀ℝ n
-        f = \ x_0 -> pickEqs $ eqs x_0 `monIndex` zeroMonomial
-
-        f'_x :: Vec n ( 𝕀ℝ n )
-        f'_x = fmap ( \ i -> pickEqs $ eqs x `monIndex` linearMonomial i ) ( universe @n )
-
-  = let -- Interval Newton method: take one Gauss–Seidel step
-        -- for the system of equations f'(x) ( x - x_mid ) = - f(x_mid).
-        minus_f_x_mid = unT $ -1 *^ T ( boxMidpoint $ f x_mid )
-
-        -- Precondition the above linear system into A ( x - x_mid ) = B.
-        ( a, b ) = precondition precondMeth
-                      f'_x ( singleton minus_f_x_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 ( \ x' -> unT $ x' ^+^ T x_mid ) )
-                  $ gaussSeidelUpdate gsUpdate a b ( T x ^-^ T x_mid )
-        ( done, todo ) = bimap ( map fst ) ( map fst )
-                       $ partition snd gsGuesses
-    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.
-      do tell $ noDoneBoxes { doneSolBoxes = done }
-         return todo
-{-# INLINEABLE intervalGaussSeidel #-}
-{-
-
-mbDeg = topologicalDegree 0.005 f x
-det = case f'_x of
-        Vec [ c1, c2 ] ->
-          let a_11 = c1 `index` Fin 1
-              a_12 = c2 `index` Fin 1
-              a_21 = c1 `index` Fin 2
-              a_22 = c2 `index` Fin 2
-          in a_11 * a_22 - a_12 * a_21
-        _ -> error "TODO: just testing n=2 here"
-
-if | not $ 0 ∈ det
-   , mbDeg == Just 0
-   -> return []
-     -- If the Jacobian is invertible over the box, then the topological
-     -- degree tells us exactly how many solutions there are in the box.
--}
-
 -- | A partial or complete Gauss–Seidel step for the equation \( A X = B \),
 -- refining the initial guess box for \( X \) into up to \( 2^n \) (disjoint) new boxes.
 gaussSeidelUpdate
@@ -291,22 +190,6 @@ gaussSeidelStep_Complete as b ( T x0 ) = coerce $ do
   return ( x', and subs )
 {-# INLINEABLE gaussSeidelStep_Complete #-}
 
-fromComponents
-  :: forall n
-  .  ( Representable Double ( ℝ n ), n ~ RepDim ( ℝ n ) )
-  => ( Fin n -> [ ( 𝕀 Double, Bool ) ] ) -> [ ( 𝕀ℝ n, Vec n Bool ) ]
-fromComponents f = do
-  ( xs, bs ) <- unzip <$> traverse f ( universe @n )
-  return $ ( tabulate $ \ i -> xs ! i, bs )
-    -- TODO: this could be more efficient.
-{-# INLINEABLE fromComponents #-}
-
--- | The midpoint of a box.
-boxMidpoint :: Representable Double ( ℝ n ) => 𝕀ℝ n -> ℝ n
-boxMidpoint box =
-  tabulate $ \ i -> midpoint ( box `index` i )
-{-# INLINEABLE boxMidpoint #-}
-
 -- | Pre-condition the system \( AX = B \).
 precondition
   :: forall n
@@ -340,16 +223,3 @@ precondition meth as b =
         )
         ( 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/Newton/LP.hs

@@ -0,0 +1,132 @@
+{-# LANGUAGE CApiFFI #-}
+{-# LANGUAGE ForeignFunctionInterface #-}
+{-# LANGUAGE ParallelListComp #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+-- | A linear programming approach for solving systems of
+-- interval linear equations.
+module Math.Root.Isolation.Newton.LP
+  ( solveIntervalLinearEquations )
+  where
+
+-- base
+import Control.Arrow
+  ( (&&&) )
+import Data.Coerce
+  ( coerce )
+import Data.Foldable
+  ( toList )
+import Foreign.C.Types
+  ( CDouble(..), CInt(..), CUInt(..) )
+import Foreign.Marshal
+  ( allocaArray, peekArray, pokeArray, with, withArray )
+import Foreign.Ptr
+  ( Ptr, castPtr )
+import Foreign.Storable
+  ( Storable(..) )
+import GHC.TypeNats
+  ( KnownNat, natVal' )
+import GHC.Exts
+  ( proxy# )
+import System.IO.Unsafe
+  ( unsafePerformIO )
+
+-- MetaBrush
+import Math.Interval
+import Math.Linear
+
+--------------------------------------------------------------------------------
+-- Linear programming approach to solving systems of linear equations
+-- (in two variables).
+
+-- | Solve the system of linear equations \( A X = B \)
+-- using linear programming.
+solveIntervalLinearEquations
+  :: ( KnownNat d, Representable Double ( ℝ d ), Show ( ℝ d ) )
+  => Vec 2 ( 𝕀ℝ d )    -- ^ columns of \( A \)
+  -> 𝕀ℝ d              -- ^ \( B \)
+  -> T ( 𝕀ℝ 2 )        -- ^ initial box \( X \)
+  -> [ ( T ( 𝕀ℝ 2 ), Bool ) ]
+solveIntervalLinearEquations a b x =
+  let !sols = unsafePerformIO $ intervalSystem2D_LP a b x
+  in if
+      | any hasNaNs sols
+      -> [ ( x, False ) ]
+      | length sols <= 1
+      -> map ( id &&& isSubBox x ) sols
+      | otherwise
+      -> map ( , False ) sols
+
+-- Assuming the second box is a subset of the second box, returns whether it
+-- is in fact a strict subset.
+isSubBox :: T ( 𝕀ℝ 2 ) -> T ( 𝕀ℝ 2 ) -> Bool
+isSubBox
+  ( T ( 𝕀 ( ℝ2 x1_min y1_min ) ( ℝ2 x1_max y1_max ) ) )
+  ( T ( 𝕀 ( ℝ2 x2_min y2_min ) ( ℝ2 x2_max y2_max ) ) )
+    =  ( ( x2_min > x1_min && x2_max < x1_max ) || x2_min == x2_max )
+    && ( ( y2_min > y1_min && y2_max < y1_max ) || y2_min == y2_max )
+
+hasNaNs :: T (𝕀ℝ 2) -> Bool
+hasNaNs ( T ( 𝕀 ( ℝ2 x_min y_min ) ( ℝ2 x_max y_max ) ) ) =
+  any ( \ x -> isNaN x || isInfinite x ) [ x_min, y_min, x_max, y_max ]
+
+intervalSystem2D_LP
+  :: forall d
+  .  ( KnownNat d, Representable Double ( ℝ d ) )
+  => Vec 2 ( 𝕀ℝ d )
+  -> 𝕀ℝ d
+  -> T ( 𝕀ℝ 2 )
+  -> IO [ T ( 𝕀ℝ 2 ) ]
+intervalSystem2D_LP a b x =
+  allocaArray 4 \ ptrSolutions ->
+    with ( CBox $ unT $ x ) \ ptrBox ->
+      withArray ( mkEquationArray a b ) \ ptrEqs -> do
+        CInt nbSols <-
+          interval_system_2d ptrSolutions ptrBox ptrEqs ( fromIntegral d )
+        if nbSols < 0
+        then
+          error $ unlines
+            [ "interval_system_2d returned with exit code " ++ show nbSols
+            , "This probably means it was given invalid input." ]
+        else
+          coerce <$> peekArray ( fromIntegral nbSols ) ptrSolutions
+  where
+    d = natVal' @d proxy#
+
+mkEquationArray
+  :: ( KnownNat d, Representable Double ( ℝ d ) )
+  => Vec 2 ( 𝕀ℝ d ) -> 𝕀ℝ d -> [ CEqn ]
+mkEquationArray ( Vec [ a_x, a_y ] ) b =
+  [ CEqn a_x_i a_y_i b_i
+  | a_x_i <- toList $ coordinates a_x
+  | a_y_i <- toList $ coordinates a_y
+  | b_i <- toList $ coordinates b
+  ]
+mkEquationArray _ _ = error "impossible"
+
+foreign import ccall "interval_system_2d"
+  interval_system_2d :: Ptr CBox -> Ptr CBox -> Ptr CEqn -> CUInt -> IO CInt
+
+data CEqn = CEqn !( 𝕀 Double ) !( 𝕀 Double ) !( 𝕀 Double )
+instance Storable CEqn where
+  sizeOf _ = 6 * sizeOf @Double undefined
+  alignment _ = 4 * alignment @Double undefined
+  peek ptr = do
+    [ CDouble a_min, CDouble a_max, CDouble b_min, CDouble b_max, CDouble c_min, CDouble c_max ]
+      <- peekArray 6 ( castPtr ptr :: Ptr CDouble )
+    return $
+      CEqn ( 𝕀 a_min a_max ) ( 𝕀 b_min b_max ) ( 𝕀 c_min c_max )
+  poke ptr ( CEqn ( 𝕀 a_min a_max ) ( 𝕀 b_min b_max ) ( 𝕀 c_min c_max ) )
+    = pokeArray ( castPtr ptr ) [ CDouble a_min, CDouble a_max, CDouble b_min, CDouble b_max, CDouble c_min, CDouble c_max ]
+
+
+newtype CBox = CBox ( 𝕀ℝ 2 )
+instance Storable CBox where
+  sizeOf _ = 4 * sizeOf @Double undefined
+  alignment _ = 4 * alignment @Double undefined
+  peek ptr = do
+    [ CDouble x_min, CDouble x_max, CDouble y_min, CDouble y_max ] <- peekArray 4 ( castPtr ptr :: Ptr CDouble )
+    return $
+      CBox ( 𝕀 ( ℝ2 x_min y_min ) ( ℝ2 x_max y_max ) )
+  poke ptr ( CBox ( 𝕀 ( ℝ2 x_min y_min ) ( ℝ2 x_max y_max ) ) ) =
+    pokeArray ( castPtr ptr ) [ CDouble x_min, CDouble x_max, CDouble y_min, CDouble y_max ]

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

@@ -0,0 +1,48 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+
+-- | Utilities for root isolation
+module Math.Root.Isolation.Utils
+  ( fromComponents, matMulVec, boxMidpoint )
+  where
+
+-- base
+import Prelude hiding ( unzip )
+import Data.Foldable
+  ( toList )
+import Data.List.NonEmpty
+  ( unzip )
+
+-- MetaBrush
+import Math.Interval
+import Math.Linear
+
+--------------------------------------------------------------------------------
+
+fromComponents
+  :: forall n
+  .  ( Representable Double ( ℝ n ), n ~ RepDim ( ℝ n ) )
+  => ( Fin n -> [ ( 𝕀 Double, Bool ) ] ) -> [ ( 𝕀ℝ n, Vec n Bool ) ]
+fromComponents f = do
+  ( xs, bs ) <- unzip <$> traverse f ( universe @n )
+  return $ ( tabulate $ \ i -> xs ! i, bs )
+    -- TODO: this could be more efficient.
+{-# INLINEABLE fromComponents #-}
+
+-- | The midpoint of a box.
+boxMidpoint :: Representable Double ( ℝ n ) => 𝕀ℝ n -> ℝ n
+boxMidpoint box =
+  tabulate $ \ i -> midpoint ( box `index` i )
+{-# INLINEABLE boxMidpoint #-}
+
+-- | 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 #-}