Fix extendedRecip (negative infinity)

brush-strokes/bench/Main.hs

@@ -17,16 +17,22 @@ module Main
   where
 
 -- base
+import Control.Concurrent.MVar
+  ( newMVar )
 import Data.Coerce
   ( coerce )
 import Data.Foldable
   ( for_ )
+import Data.List
+  ( intercalate )
 import GHC.Exts
   ( Proxy#, proxy# )
 import GHC.Generics
   ( Generic )
 import GHC.TypeNats
   ( type (-) )
+import Numeric
+  ( showFFloat )
 
 -- containers
 import Data.Sequence
@@ -38,6 +44,8 @@ import Data.Tree
 
 -- brush-strokes
 import Calligraphy.Brushes
+import Debug.Utils
+  ( logToFile )
 import Math.Algebra.Dual
 import Math.Bezier.Spline
 import Math.Bezier.Stroke
@@ -66,6 +74,9 @@ main = for_ testCases $ \ testCase@( TestCase { testName, testAlgorithmParams }
       , " dunno: " ++ show dunno
       , "  sols: " ++ show sols
       ]
+  --logFileMVar <- newMVar "logs/trickyCusp.log"
+  --logToFile logFileMVar (unlines logLines)
+  --  `seq` return ()
 
 testCases :: [ TestCase ]
 testCases = [ ellipse , trickyCusp2 ]
@@ -140,9 +151,96 @@ putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees $ mkBox (0.5798, 0.
 ([ℝ1 0.5798, ℝ1 0.675],3,[ℝ1 0.26798, ℝ1 0.26799]) (area 0.000001) N []
  └─ ([ℝ1 0.5973000285624527, ℝ1 0.6750000000000002],3,[ℝ1 0.26798, ℝ1 0.26799000000000006]) (area 0.000001) NoSolution "ee" ([ℝ1 0.5973000285624527, ℝ1 0.6750000000000002],3,[ℝ1 0.26798, ℝ1 0.26799000000000006])
 
+eval fI $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799)
+> D12 {_D12_v = T[ℝ2 -10088.6674944889 -3281.3820867312834, ℝ2 4124.668381545453 4524.807156085763], _D12_dx = TT[ℝ2 -173746.97965005718 -33281.18494907289, ℝ2 298.2609121556852 23639.772884799597], _D12_dy = TT[ℝ2 -18454.27716258352 -28337.509817580823, ℝ2 1163.6949532017436 -13936.383137525536]}}
+i.e.
+> f = [ℝ2 -10088.6674944889 -3281.3820867312834, ℝ2 4124.668381545453 4524.807156085763]
+> f_t = [ℝ2 -173746.97965005718 -33281.18494907289, ℝ2 298.2609121556852 23639.772884799597]
+> f_s = [ℝ2 -18454.27716258352 -28337.509817580823, ℝ2 1163.6949532017436 -13936.383137525536]
+
+(f, fI) = testCaseStrokeFunctions trickyCusp2
+t = 𝕀 (ℝ1 0.5798) (ℝ1 0.675)
+s = 𝕀 (ℝ1 0.26798) (ℝ1 0.26799)
+t_mid = 0.5 * ( 0.5798 + 0.675 )
+s_mid = 0.5 * ( 0.26798 + 0.26799 )
+D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) = dEdsdcdt $ eval fI (t, 3, s)
+t' = coerce ( (-) @( 𝕀 Double ) ) t ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1
+s' = coerce ( (-) @( 𝕀 Double ) ) s ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1
+a = ( f_t, f_s )
+b = negV2 $ singleton $ midV2 f
+[((t2', s2'), isContr)] = gaussSeidel a b (t', s')
+t2 = coerce ( (+) @( 𝕀 Double ) ) t2' ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1
+s2 = coerce ( (+) @( 𝕀 Double ) ) s2' ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1
+
+t2
+> [ℝ1 0.6102365832093095, ℝ1 0.6750000000000002]
+s2
+> [ℝ1 0.26798, ℝ1 0.26799000000000006]
+
+
+let ( 𝕀 ( ℝ2 a11_lo a21_lo ) ( ℝ2 a11_hi a21_hi ), 𝕀 ( ℝ2 a12_lo a22_lo ) ( ℝ2 a12_hi a22_hi ) ) = a
+let ( 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi ) ) = b
+let ( 𝕀 ( ℝ1 x1_lo ) ( ℝ1 x1_hi ), 𝕀 ( ℝ1 x2_lo ) ( ℝ1 x2_hi ) ) = ( t', s' )
+
+a11 = 𝕀 a11_lo a11_hi
+a12 = 𝕀 a12_lo a12_hi
+a21 = 𝕀 a21_lo a21_hi
+a22 = 𝕀 a22_lo a22_hi
+b1  = 𝕀  b1_lo  b1_hi
+b2  = 𝕀  b2_lo  b2_hi
+x1  = 𝕀  x1_lo  x1_hi
+x2  = 𝕀  x2_lo  x2_hi
+
+( b1 - a12 * x2 )
+> [2981.90728508591, 2982.0918278575364]
+
+extendedRecip a11
 
 -}
 
+negV2 :: 𝕀ℝ 2 -> 𝕀ℝ 2
+negV2 ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) =
+  let !( 𝕀 x'_lo x'_hi ) = negate $ 𝕀 x_lo x_hi
+      !( 𝕀 y'_lo y'_hi ) = negate $ 𝕀 y_lo y_hi
+  in 𝕀 ( ℝ2 x'_lo y'_lo ) ( ℝ2 x'_hi y'_hi )
+
+midV2 :: 𝕀ℝ 2 -> ℝ 2
+midV2 ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) =
+  ℝ2 ( 0.5 * ( x_lo + x_hi ) ) ( 0.5 * ( y_lo + y_hi ) )
+
+logLines :: [ String ]
+logLines =
+  [ "f = dE/ds * dc/dt: f, df/dt, df/ds"
+  , "{" ++
+    (intercalate ","
+     [ "{" ++ showD (midPoint t) ++ "," ++ showD (midPoint s) ++ ",{" ++ intercalate "," vals ++ "}}"
+     | t <- map ( around 0.5798  ) [-0.05, -0.049.. 0.05]
+     , let i = 3
+     , s <- map ( around 0.26798 ) [-0.05, -0.049.. 0.05]
+     , let StrokeDatum
+             { 𝛿E𝛿sdcdt = D12 (T f) (T (T f_t)) (T (T f_s))
+             } = (curvesI t `Seq.index` i) s
+           ℝ2 vx vy = midPoint2 f
+           ℝ2 vx_t vy_t = midPoint2 f_t
+           ℝ2 vx_s vy_s = midPoint2 f_s
+           vals = [ "{" ++ showD vx ++ "," ++ showD vy ++ "}"
+                  , "{" ++ showD vx_t ++ "," ++ showD vy_t ++ "}"
+                  , "{" ++ showD vx_s ++ "," ++ showD vy_s ++ "}"
+                  ]
+     ]
+    ) ++ "}"
+  ]
+  where
+    around :: Double -> Double -> 𝕀ℝ 1
+    around z0 z = 𝕀 ( ℝ1 ( z + z0 - 1e-6 ) ) ( ℝ1 ( z + z0 + 1e-6 ) )
+    ( _, curvesI ) = testCaseStrokeFunctions trickyCusp2
+    midPoint (𝕀 (ℝ1 lo) (ℝ1 hi)) = 0.5 * ( lo + hi )
+    midPoint2 (𝕀 (ℝ2 lo_x lo_y) (ℝ2 hi_x hi_y))
+      = ℝ2 ( 0.5 * ( lo_x + hi_x ) ) ( 0.5 * ( lo_y + hi_y ) )
+
+showD :: Double -> String
+showD float = showFFloat (Just 6) float ""
+
 --------------------------------------------------------------------------------
 
 ellipse :: TestCase

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

@@ -24,7 +24,7 @@ module Math.Bezier.Stroke
   , IntervalNewtonStep(..)
   , IntervalNewtonLeaf(..)
   , Box
-  , intervalNewtonGS, intervalNewtonGSFrom
+  , intervalNewtonGS, intervalNewtonGSFrom, gaussSeidel
 
   )
   where
@@ -1288,18 +1288,6 @@ intersect ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 )
     lo = max lo1 lo2
     hi = min hi1 hi2
 
-extendedDivide :: 𝕀 Double -> 𝕀 Double -> [ 𝕀 Double ]
-extendedDivide x y = map ( x * ) ( extendedRecip y )
-
-extendedRecip :: 𝕀 Double -> [ 𝕀 Double ]
-extendedRecip x@( 𝕀 lo hi )
-  | lo == 0 && hi == 0
-  = [ 𝕀 ( -1 / 0 ) ( 1 / 0 ) ]
-  | lo >= 0 || hi <= 0
-  = [ recip x ]
-  | otherwise
-  = [ recip ( 𝕀 lo 0 ), recip ( 𝕀 0 hi ) ]
-
 -- | Computes the brush stroke coordinates of a cusp from
 -- the @(t,s)@ parameter values.
 cuspCoords :: ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) )
@@ -1421,7 +1409,7 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox =
         <- ( 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 ) ) }
+                    , 𝛿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
       , let ee_potential_zero = inf ee <= ℝ1 0 && sup ee >= ℝ1 0
             𝛿E𝛿sdcdt_potential_zero = cmpℝ2 (<=) ( inf f ) ( ℝ2 0 0 ) && cmpℝ2 (>=) ( sup f ) ( ℝ2 0 0 )
@@ -1433,7 +1421,7 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox =
                   -- for the equation f'(x) v = - f(x_mid),
                   -- where f = 𝛿E/𝛿s * dc/dt
                   !( a, b ) = precondition precondMethod
-                                ( f_t_mid, f_s_mid )
+                                ( midI f_t, midI f_s )
                                 ( f_t, f_s ) ( neg f_mid )
                   --(a, b)
                   --  | 𝕀 (ℝ1 ee_lo) (ℝ1 ee_hi) <- ee_mid
@@ -1482,6 +1470,9 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox =
            | otherwise
            -> return [ IntervalNewtonLeaf $ NoSolution (if ee_potential_zero then "dc/dt" else "ee") cand ]
       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 ) )
         t_mid = 0.5 * ( t_lo + t_hi )
         s_mid = 0.5 * ( s_lo + s_hi )
         i_t_mid = singleton ( ℝ1 t_mid )

brush-strokes/src/Math/Interval.hs

@@ -11,11 +11,12 @@ module Math.Interval
   , singleton, nonDecreasing
   , inside
   , aabb
+  , extendedDivide, extendedRecip
   )
   where
 
 -- base
-import Prelude hiding ( Num(..) )
+import Prelude hiding ( Num(..), Fractional(..) )
 
 -- acts
 import Data.Act
@@ -69,6 +70,26 @@ instance Act   ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where
 instance Torsor ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where
   a --> b = T $ b - a
 
+-------------------------------------------------------------------------------
+-- Extended division
+
+extendedDivide :: ( Field d, Field ( 𝕀 d ), Ord d ) => 𝕀 d -> 𝕀 d -> [ 𝕀 d ]
+extendedDivide x y = map ( x * ) ( extendedRecip y )
+{-# SPECIALISE extendedDivide :: 𝕀 Double -> 𝕀 Double -> [ 𝕀 Double ] #-}
+
+extendedRecip :: ( Field d, Field ( 𝕀 d ), Ord d ) => 𝕀 d -> [ 𝕀 d ]
+extendedRecip x@( 𝕀 lo hi )
+  | lo == fromInteger 0 && hi == fromInteger 0
+  = [ 𝕀 negInf posInf ]
+  | lo >= fromInteger 0 || hi <= fromInteger 0
+  = [ recip x ]
+  | otherwise
+  = [ 𝕀 negInf ( recip lo ), 𝕀 ( recip hi ) posInf ]
+  where
+    negInf = fromInteger (-1) / fromInteger 0
+    posInf = fromInteger   1  / fromInteger 0
+{-# SPECIALISE extendedRecip :: 𝕀 Double -> [ 𝕀 Double ] #-}
+
 -------------------------------------------------------------------------------
 -- Lattices.
 

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

@@ -109,7 +109,11 @@ instance Prelude.Fractional ( 𝕀 Double ) where
     in 𝕀 q q
   recip (𝕀 lo hi)
 -- #ifdef ASSERTS
-    | lo >= 0 || hi <= 0
+    | lo == 0
+    = 𝕀 ( fst $ divI 1 hi ) ( 1 / 0 )
+    | hi == 0
+    = 𝕀 ( -1 / 0 ) ( snd $ divI 1 lo )
+    | lo > 0 || hi < 0
 -- #endif
     = 𝕀 ( fst $ divI 1 hi ) ( snd $ divI 1 lo )
 -- #ifdef ASSERTS
@@ -118,7 +122,11 @@ instance Prelude.Fractional ( 𝕀 Double ) where
 -- #endif
   𝕀 x_lo x_hi / 𝕀 y_lo y_hi
 -- #ifdef ASSERTS
-    | y_lo >= 0 || y_hi <= 0
+    | y_lo == 0
+    = 𝕀 ( fst $ divI x_lo y_hi ) ( 1 / 0 )
+    | y_hi == 0
+    = 𝕀 ( -1 / 0 ) ( snd $ divI x_hi y_lo )
+    | y_lo > 0 || y_hi < 0
 -- #endif
     = 𝕀 ( fst $ divI x_lo y_hi ) ( snd $ divI x_hi y_lo )
 -- #ifdef ASSERTS