Implement box(1) consistency check for cusp finding

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

@@ -610,32 +610,13 @@ x2  = 𝕀  x2_lo  x2_hi
 -}
 
 {-
-SMOKING GUN:
 
-> (f, fI) = brushStrokeFunctions $ ellipseBrushStroke (0,1) pi
-> (t, i, s) = mkBox (0.54, 0.55) 2 (0.5480, 0.5481)
-> _D32_dxdxdx $ dbrush $ eval fI (t,i,s)
-> [ℝ2 54187.61174031432 -11626.47655724034, ℝ2 55257.54384135226 -9759.09092706883]
+bound consistency / box consistency
 
-In Mathematica:
+ghci> (f, fI) = brushStrokeFunctions $ ellipseBrushStroke (0,1) pi
+ghci> ee $ eval fI $ mkBox (0, succFP 0) 2 (0.5,0.75)
+D22 {_D22_v = [ℝ1 -4379.198802308688, ℝ1 -1791.4548164846506], _D22_dx = T[ℝ1 -3417.0023260831304, ℝ1 6215.362085329192], _D22_dy = T[ℝ1 3524.130928930924, ℝ1 8849.881134814816], _D22_dxdx = T[ℝ1 -13397.019488083239, ℝ1 61439.0409103
 
-  FindMinimum[{D[b[vt, s, myPb0, myPb1][[1]], {vt, 3}] /. (vt -> t), 
-    0.54 <= t <= 0.55 && 0.548 <= s <= 0.5481}, {{t, 0.545}, {s, 
-     0.54805}}]
-> { -503.445, {t -> 0.54, s -> 0.548} }
 
-  FindMaximum[{D[b[vt, s, myPb0, myPb1][[1]], {vt, 3}] /. (vt -> t), 
-    0.54 <= t <= 0.55 && 0.548 <= s <= 0.5481}, {{t, 0.545}, {s, 
-     0.54805}}]
-> { -480.086, {t -> 0.55, s -> 0.5481} }
-
-so the actual range is [-503.445,-480.086] but we have computed a range of
-[54187.61174031432, 55257.54384135226] which does not contain the actual range!
-
--}
-
-{-
-
-Math.Ring.cos  ( linearD @_ @3 (\i -> 𝕀 (getR1 $ inf i) (getR1 $ sup i)) (𝕀 (ℝ1 $ 0.548 * pi) (ℝ1 $ 0.5481 * pi)) :: D3𝔸1 ( 𝕀 Double ) )
 
 -}

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

@@ -35,7 +35,7 @@ import Control.Arrow
 import Control.Concurrent.MVar
   ( MVar, newMVar )
 import Control.Monad
-  ( unless )
+  ( unless, when )
 import Control.Monad.ST
   ( RealWorld, ST )
 import Data.Bifunctor
@@ -58,6 +58,8 @@ import Data.Maybe
   ( fromMaybe, isJust, listToMaybe, mapMaybe )
 import Data.Semigroup
   ( sconcat )
+import Data.Traversable
+  ( for )
 import GHC.Exts
   ( newMutVar#, runRW#, inline
   , Proxy#, proxy#
@@ -108,8 +110,8 @@ import Control.Monad.Trans.Class
   ( lift )
 import Control.Monad.Trans.Except
   ( Except, runExcept, throwE )
-import Control.Monad.Trans.State.Strict
-  ( StateT, runStateT, evalStateT, get, put )
+import Control.Monad.Trans.State.Strict as State
+  ( StateT, State, runStateT, evalState, evalStateT, get, put )
 import Control.Monad.Trans.Writer.CPS
   ( Writer, WriterT, execWriterT, runWriter, tell )
 
@@ -138,6 +140,8 @@ import Math.Differentiable
   ( Differentiable, DiffInterp, I )
 import Math.Epsilon
   ( epsilon, nearZero )
+import Math.Float.Utils
+  ( succFP, prevFP )
 import Math.Interval
 import Math.Linear
 import Math.Module
@@ -1493,7 +1497,7 @@ intervalNewtonGSFrom
   -> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) )
 intervalNewtonGSFrom precondMethod minWidth eqs initBox =
   runWriter $
-    map ( initBox , ) <$> evalStrokeDataAndGo initBox
+    map ( initBox , ) <$> go initBox
 
   where
 
@@ -1505,19 +1509,13 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox =
       rest <- traverse ( \ c -> do { trees <- r c; return [ (c, tree) | tree <- trees ] } ) cands
       return [ f $ concat rest ]
 
-    evalStrokeDataAndGo :: Box -> Writer ( [Box], [Box] ) [ IntervalNewtonTree Box ]
-    evalStrokeDataAndGo box@( t, i, s ) = go ( box, ( eqs t `Seq.index` i ) s )
-
-    go :: ( Box, StrokeDatum 3 𝕀  ) -- box to work on
+    go :: Box -- box to work on
        -> Writer ( [ Box ], [ Box ] )
                  [ IntervalNewtonTree Box ]
-    go ( cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
-              , i
-              , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
-              )
-       , sd@( StrokeDatum { ee = D22 _ ( T ee_t ) (T ee_s ) _ _ _
-                          , 𝛿E𝛿sdcdt = D12 _ ( T ( T f_t ) ) ( T ( T f_s ) ) })
-       )
+    go cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
+            , i
+            , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
+            )
       -- Box is small: stop processing it.
       | t_hi - t_lo < minWidth && s_hi - s_lo < minWidth
       = do let res = TooSmall cand
@@ -1527,6 +1525,24 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox =
         ee_potential_zero sd
         -- Check the 𝛿E𝛿sdcdt interval contains zero.
       , 𝛿E𝛿sdcdt_potential_zero sd
+      = continue cand
+      -- Box doesn't contain a solution: discard it.
+      | otherwise
+      = return [ IntervalNewtonLeaf $ NoSolution ( if ee_potential_zero sd then "dc/dt" else "ee" ) cand ]
+      where
+        sd = ( eqs t `Seq.index` i ) s
+
+    continue :: Box
+             -> Writer ( [ Box ], [ Box ] )
+                       [ IntervalNewtonTree Box ]
+    continue
+      cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
+           , i
+           , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
+           )
+      | StrokeDatum { ee = D22 _ ( T ee_t ) (T ee_s ) _ _ _
+                    , 𝛿E𝛿sdcdt = D12 _ ( T ( T f_t ) ) ( T ( T f_s ) ) }
+        <- ( eqs t `Seq.index` i ) s
       , StrokeDatum { ee = D22 _ee_mid _ _ _ _ _
                     , 𝛿E𝛿sdcdt = D12 ( T f_mid ) ( T ( T _f_t_mid ) ) ( T ( T _f_s_mid ) ) }
         <- ( eqs i_t_mid `Seq.index` i ) i_s_mid
@@ -1564,7 +1580,7 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox =
              in do tell ([], done)
                    case todo of
                      [] -> return [ IntervalNewtonLeaf $ NoSolution "GaussSeidel" cand ]
-                     _ -> recur evalStrokeDataAndGo ( IntervalNewtonStep ( IntervalNewtonContraction done ) )
+                     _ -> recur go ( IntervalNewtonStep ( IntervalNewtonContraction done ) )
                             todo
            else
               -- Gauss–Seidel failed to shrink the boxes, so bisect instead.
@@ -1605,40 +1621,23 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox =
                    | d_skip && u_skip
                    = ( [], ( "du", s_mid ) )
                    | l_skip
-                   = ( [ ( r, r_dat ) ], ( "r", t_mid ) )
+                   = ( [ r ], ( "r", t_mid ) )
                    | r_skip
-                   = ( [ ( l, l_dat ) ], ( "l", t_mid ) )
+                   = ( [ l ], ( "l", t_mid ) )
                    | d_skip
-                   = ( [ ( u, u_dat ) ], ( "u", s_mid ) )
+                   = ( [ u ], ( "u", s_mid ) )
                    | u_skip
-                   = ( [ ( d, d_dat ) ], ( "d", s_mid ) )
+                   = ( [ d ], ( "d", s_mid ) )
                    | tJWidth >= sJWidth
-                 --  t_hi - t_lo >= s_hi - s_lo
-                   = ( [ ( l, l_dat ), ( r, r_dat ) ], ( "t", t_mid ) )
+                   -- ... but don't allow thin boxes.
+                   || ( t_hi - t_lo >= 10 * ( s_hi - s_lo ) )
+                   && not ( s_hi - s_lo >= 10 * ( t_hi - t_lo ) )
+                   = ( [ l, r ], ( "t", t_mid ) )
                    | otherwise
-                   = ( [ ( d, d_dat ), ( u, u_dat ) ], ( "s", s_mid ) )
-             in recur go ( IntervalNewtonStep ( IntervalNewtonBisection whatBis ) . map (first fst) )
-                  bisGuesses
-
-           -- Box doesn't contain a solution: discard it.
-           | otherwise
-           = return [ IntervalNewtonLeaf $ NoSolution ( if ee_potential_zero sd then "dc/dt" else "ee" ) cand ]
+                   = ( [ d, u ], ( "s", s_mid ) )
+             in recur go ( IntervalNewtonStep ( IntervalNewtonBisection whatBis ) )
+                  ( makeBox1Consistent =<< bisGuesses )
       where
-        midI :: 𝕀ℝ 2 -> 𝕀ℝ 2
-        midI ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) =
-          singleton $ ℝ2 ( 0.5 * ( x_lo + x_hi ) ) ( 0.5 * ( y_lo + y_hi ) )
-        --width :: 𝕀ℝ 1 -> Double
-        --width ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = hi - lo
-        --normI :: 𝕀ℝ 1 -> Double
-        --normI ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = sqrt $ sup $ ( 𝕀 lo hi ) Ring.^ 2
-        --normVI :: 𝕀ℝ 2 -> Double
-        --normVI ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) =
-        --  sqrt $ sup $ ( 𝕀 x_lo x_hi ) Ring.^ 2 + ( 𝕀 y_lo y_hi ) Ring.^ 2
-        normVI3 :: 𝕀ℝ 1 -> 𝕀ℝ 2 -> Double
-        normVI3 ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) )
-          = sqrt $ max ( abs lo   ) ( abs hi   ) Ring.^ 2
-                 + max ( abs x_lo ) ( abs x_hi ) Ring.^ 2
-                 + max ( abs y_lo ) ( abs y_hi ) Ring.^ 2
         t_mid = 0.5 * ( t_lo + t_hi )
         s_mid = 0.5 * ( s_lo + s_hi )
         i_t_mid = singleton ( ℝ1 t_mid )
@@ -1656,6 +1655,19 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox =
                 !( 𝕀 y'_lo y'_hi ) = negate $ 𝕀 y_lo y_hi
             in 𝕀 ( ℝ2 x'_lo y'_lo ) ( ℝ2 x'_hi y'_hi )
 
+    -- 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 :: Box -> [ Box ]
+    makeBox1Consistent x =
+      ( `State.evalState` False ) $
+        pipeFunctionsWhileTrue (allNarrowingOperators eps_eq eps_bis eqs) x
+    eps_eq = 1e-5 -- TODO: this is an additive constant,
+                  -- but I think we also want a multiplicative constant?
+    eps_bis = minWidth
+
 ee_potential_zero :: StrokeDatum 3 𝕀 -> Bool
 ee_potential_zero dat =
      inf ( _D22_v $ ee dat ) <= ℝ1 0
@@ -1716,3 +1728,172 @@ matMulVec ( ℝ2 a11 a21, ℝ2 a12 a22 ) ( 𝕀 ( ℝ2 u_lo v_lo ) ( ℝ2 u_hi v
 cmpℝ2 :: ( Double -> Double -> Bool ) -> ℝ 2 -> ℝ 2 -> Bool
 cmpℝ2 cmp ( ℝ2 x1 y1 ) ( ℝ2 x2 y2 )
   = cmp x1 x2 && cmp y1 y2
+
+midI :: 𝕀ℝ 2 -> 𝕀ℝ 2
+midI ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) =
+  singleton $ ℝ2 ( 0.5 * ( x_lo + x_hi ) ) ( 0.5 * ( y_lo + y_hi ) )
+--width :: 𝕀ℝ 1 -> Double
+--width ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = hi - lo
+--normI :: 𝕀ℝ 1 -> Double
+--normI ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = sqrt $ sup $ ( 𝕀 lo hi ) Ring.^ 2
+--normVI :: 𝕀ℝ 2 -> Double
+--normVI ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) =
+--  sqrt $ sup $ ( 𝕀 x_lo x_hi ) Ring.^ 2 + ( 𝕀 y_lo y_hi ) Ring.^ 2
+normVI3 :: 𝕀ℝ 1 -> 𝕀ℝ 2 -> Double
+normVI3 ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) =
+  sqrt $ max ( abs lo   ) ( abs hi   ) Ring.^ 2
+       + max ( abs x_lo ) ( abs x_hi ) Ring.^ 2
+       + max ( abs y_lo ) ( abs y_hi ) Ring.^ 2
+
+-- 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 eps_equal eps_bisection 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 delta <- f_x_left `extendedDivide` f'_x --f'_x_left
+                       let x_new = x_left - delta
+                       map fst $ x_new `intersect` x
+          in
+            if | null x's
+               -> []
+               | ( width x - sum ( map width x's ) ) < eps_equal
+               -- TODO: do a check with the relative reduction in size?
+               -> x's
+               | otherwise
+               -> do
+                x' <- x's
+                if sup x' - inf x' < eps_bisection
+                then return x'
+                else left_narrow =<< bisect x'
+
+-- TODO: de-duplicate with 'leftNarrow'?
+rightNarrow :: Double
+            -> Double
+            -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) )
+            -> 𝕀 Double
+            -> [ 𝕀 Double ]
+rightNarrow eps_equal eps_bisection 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 delta <- f_x_right `extendedDivide` f'_x --f'_x_right
+                       let x_new = x_right - delta
+                       map fst $ x_new `intersect` x
+          in
+            if | null x's
+               -> []
+               | ( width x - sum ( map width x's ) ) < eps_equal
+               -- TODO: do a check with the relative reduction in size?
+               -> x's
+               | otherwise
+               -> do
+                x' <- x's
+                if sup x' - inf x' < eps_bisection
+                then return x'
+                else right_narrow =<< bisect x'
+
+bisect :: 𝕀 Double -> [ 𝕀 Double ]
+bisect x@( 𝕀 x_inf x_sup )
+  | x_inf == x_sup
+  = [ x ]
+  | otherwise
+  = [ 𝕀 x_inf x_mid, 𝕀 x_mid x_sup ]
+  where x_mid = 0.5 * ( x_inf + x_sup )
+
+-- | 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
+
+allNarrowingOperators
+  :: Double
+  -> Double
+  -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
+  -> [ Box -> State Bool [ Box ] ]
+allNarrowingOperators eps_eq eps_bis eqs =
+  [ \ cand ->
+    let newCands = map ( `setter` cand ) $ narrowFn eps_eq eps_bis ( \ x0 -> ff' $ setter x0 cand ) ( getter cand )
+    in do
+      when ( ( width ( getter cand ) - sum ( map ( width . getter ) newCands ) ) > eps_eq ) $
+        -- TODO: do a check with the relative reduction in size?
+        State.put True
+      return newCands
+  | narrowFn <- [ leftNarrow, rightNarrow ]
+  , ( ( getter, setter ), ff' ) <-
+      [ ( ( get_t, set_t ), ff' ) | ff' <- [ ff'_t, g1g1'_t, g2g2'_t ] ]
+        ++
+      [ ( ( get_s, set_s ), ff' ) | ff' <- [ ff'_s, g1g1'_s, g2g2'_s ] ]
+      ]
+  where
+    get_t ( 𝕀 ( ℝ1 t_inf ) ( ℝ1 t_sup ), _, _ ) = 𝕀 t_inf t_sup
+    set_t ( 𝕀 t_inf t_sup ) ( _, i, s ) = ( 𝕀 ( ℝ1 t_inf ) ( ℝ1 t_sup ), i, s )
+    get_s ( _ , _, 𝕀 ( ℝ1 s_inf ) ( ℝ1 s_sup ) ) = 𝕀 s_inf s_sup
+    set_s ( 𝕀 s_inf s_sup ) ( t, i, _ ) = ( t, i, 𝕀 ( ℝ1 s_inf ) ( ℝ1 s_sup ) )
+    ff'_t (t, i, s) =
+      let D22 { _D22_v = 𝕀 ( ℝ1 ee_inf ) ( ℝ1 ee_sup )
+              , _D22_dx = T ( 𝕀 ( ℝ1 ee'_inf ) ( ℝ1 ee'_sup ) )
+              } = ee $ ( eqs t `Seq.index` i ) s
+      in ( 𝕀 ee_inf ee_sup, 𝕀 ee'_inf ee'_sup )
+    ff'_s (t, i, s) =
+      let D22 { _D22_v = 𝕀 ( ℝ1 ee_inf ) ( ℝ1 ee_sup )
+              , _D22_dy = T ( 𝕀 ( ℝ1 ee'_inf ) ( ℝ1 ee'_sup ) )
+              } = ee $ ( eqs t `Seq.index` i ) s
+      in ( 𝕀 ee_inf ee_sup, 𝕀 ee'_inf ee'_sup )
+    g1g1'_t (t, i, s) =
+      let D12 { _D12_v = T ( 𝕀 ( ℝ2 x_inf _ ) ( ℝ2 x_sup _ ) )
+              , _D12_dx = T ( T ( 𝕀 ( ℝ2 x'_inf _ ) ( ℝ2 x'_sup _ ) ) )
+              } = 𝛿E𝛿sdcdt $ ( eqs t `Seq.index` i ) s
+      in ( 𝕀 x_inf x_sup, 𝕀 x'_inf x'_sup )
+    g1g1'_s (t, i, s) =
+      let D12 { _D12_v = T ( 𝕀 ( ℝ2 x_inf _ ) ( ℝ2 x_sup _ ) )
+              , _D12_dy = T ( T ( 𝕀 ( ℝ2 x'_inf _ ) ( ℝ2 x'_sup _ ) ) )
+              } = 𝛿E𝛿sdcdt $ ( eqs t `Seq.index` i ) s
+      in ( 𝕀 x_inf x_sup, 𝕀 x'_inf x'_sup )
+    g2g2'_t (t, i, s) =
+      let D12 { _D12_v = T ( 𝕀 ( ℝ2 _ y_inf ) ( ℝ2 _ y_sup ) )
+              , _D12_dx = T ( T ( 𝕀 ( ℝ2 _ y'_inf ) ( ℝ2 _ y'_sup ) ) )
+              } = 𝛿E𝛿sdcdt $ ( eqs t `Seq.index` i ) s
+      in ( 𝕀 y_inf y_sup, 𝕀 y'_inf y'_sup )
+    g2g2'_s (t, i, s) =
+      let D12 { _D12_v = T ( 𝕀 ( ℝ2 _ y_inf ) ( ℝ2 _ y_sup ) )
+              , _D12_dy = T ( T ( 𝕀 ( ℝ2 _ y'_inf ) ( ℝ2 _ y'_sup ) ) )
+              } = 𝛿E𝛿sdcdt $ ( eqs t `Seq.index` i ) s
+      in ( 𝕀 y_inf y_sup, 𝕀 y'_inf y'_sup )

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

@@ -6,6 +6,7 @@
 
 module Math.Interval
   ( 𝕀(𝕀), inf, sup
+  , width
   , scaleInterval
   , 𝕀ℝ
   , singleton, nonDecreasing
@@ -51,6 +52,9 @@ type instance D k ( 𝕀 v ) = D k v
 singleton :: a -> 𝕀 a
 singleton a = 𝕀 a a
 
+width :: AbelianGroup a => 𝕀 a -> a
+width ( 𝕀 lo hi ) = hi - lo
+
 -- | Turn a non-decreasing function into a function on intervals.
 nonDecreasing :: ( a -> b ) -> 𝕀 a -> 𝕀 b
 nonDecreasing f ( 𝕀 lo hi ) = 𝕀 ( f lo ) ( f hi )