Fix some more issues with interval recip

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

@@ -141,6 +141,9 @@ eval f ( t, i, s ) = ( f t `Seq.index` i ) s
 mkVal :: Double -> Int -> Double -> ( ℝ 1, Int, ℝ 1 )
 mkVal t i s = ( ℝ1 t, i, ℝ1 s )
 
+mkI :: ( Double, Double ) -> 𝕀 Double
+mkI ( lo, hi ) = 𝕀 lo hi
+
 mkBox :: ( Double, Double ) -> Int -> ( Double, Double ) -> Box
 mkBox ( t_min, t_max ) i ( s_min, s_max ) =
   ( 𝕀 ( ℝ1 t_min ) ( ℝ1 t_max ) , i, 𝕀 ( ℝ1 s_min ) ( ℝ1 s_max ) )
@@ -559,44 +562,34 @@ getR1 (ℝ1 u) = u
 
 (f, fI) = brushStrokeFunctions $ ellipseBrushStroke (0,1) pi
 nbPotentialSols box = let ( _newtTrees, ( dunno, sols ) ) = intervalNewtonGSFrom NoPreconditioning 1e-7 fI box in length dunno + length sols
-showTrees box = map ( uncurry showIntervalNewtonTree ) $ fst $ intervalNewtonGSFrom NoPreconditioning 1e-7 fI box
+showTrees box = putStrLn $ unlines $ map Data.Tree.View.showTree $ map ( uncurry showIntervalNewtonTree ) $ fst $ intervalNewtonGSFrom NoPreconditioning 1e-7 fI box
 
-(t, i, s) = mkBox (0.548610200176363, 0.5486102071493623) 2 (0.5480950215354709, 0.5480952)
-putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees (t,i,s)
+sol_t = 0.5486100729150693796677845183880669025324233347060776339185 :: Double
+sol_s = 0.5480950141859386853197594577293968665598143630958601978245 :: Double
+containsSol (t, _i, s) = getR1 (inf t) <= sol_t && getR1 (sup t) >= sol_t && getR1 (inf s) <= sol_s && getR1 (sup s) >= sol_s
 
+(t, i, s) = mkBox (0.54, 0.55) 2 (0.5480, 0.5481)
+containsSol (t, i, s)
+nbPotentialSols (t,i,s)
 
 t_mid = 0.5 * ( getR1 ( inf t ) + getR1 ( sup t ) )
 s_mid = 0.5 * ( getR1 ( inf s ) + getR1 ( sup s ) )
-D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) = dEdsdcdt $ eval fI (t, i, s)
+D12 ( T _f ) ( T ( T f_t ) ) ( T ( T f_s ) ) = dEdsdcdt $ eval fI (t, i, s)
+D0 ( T f_mid ) = dEdsdcdt $ eval f $ mkVal t_mid 2 s_mid
 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
+b = negV2 $ singleton f_mid
 [((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.548610200176363, ℝ1 0.5486102071493624]
+> [ℝ1 0.5451193766263323, ℝ1 0.545225929860598]
 s2
-> [ℝ1 0.5480950911334656, ℝ1 0.5480952000000001]
-
-mkBox (0.548610200176363, 0.5486102071493624) i (0.5480950911334656, 0.5480952000000001)
-
-t inf (no change)
-
-t sup (no change)
-
-s_inf:
-0.5480950215354709
-0.5480950911334656
-
-s_sup (no change)
-
-ghci> potentialCusp $ eval fI $ mkBox (0.548610200176363, 0.5486102071493623) 2 (0.54809502, 0.5480952)
-True
-ghci> potentialCusp $ eval fI $ mkBox (0.548610200176363, 0.5486102071493623) 2 (0.54809503, 0.5480952)
-False
+> [ℝ1 0.548, ℝ1 0.5481]
+containsSol (t2, i, s2)
+> False
 
 
 
@@ -613,4 +606,5 @@ b2  = 𝕀  b2_lo  b2_hi
 x1  = 𝕀  x1_lo  x1_hi
 x2  = 𝕀  x2_lo  x2_hi
 
--}
+( b1 - a12 * x2 ) `extendedDivide` a11
+-}

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

@@ -184,6 +184,8 @@ instance HasEnvelopeEquation 2 where
         !ee_t = c_tt × c_s + c_t × c_ts
         !ee_s = c_ts × c_s + c_t × c_ss
         !𝛿E𝛿sdcdt = ee_s *^ c_t ^-^ ee_t *^ c_s
+          -- TODO: we get c_t * c_t and c_s * c_s terms...
+          -- These could be squares (better with interval arithmetic)?
     in ( D12 ( co ee ) ( T $ co ee_t ) ( T $ co ee_s )
        , D0 𝛿E𝛿sdcdt )
     -- Computation of total derivative dc/dt:

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

@@ -73,22 +73,21 @@ instance Torsor ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where
 -------------------------------------------------------------------------------
 -- Extended division
 
-extendedDivide :: ( Field d, Field ( 𝕀 d ), Ord d ) => 𝕀 d -> 𝕀 d -> [ 𝕀 d ]
+extendedDivide :: 𝕀 Double -> 𝕀 Double -> [ 𝕀 Double ]
 extendedDivide x y = map ( x * ) ( extendedRecip y )
-{-# SPECIALISE extendedDivide :: 𝕀 Double -> 𝕀 Double -> [ 𝕀 Double ] #-}
 
-extendedRecip :: ( Field d, Field ( 𝕀 d ), Ord d ) => 𝕀 d -> [ 𝕀 d ]
+extendedRecip :: 𝕀 Double -> [ 𝕀 Double ]
 extendedRecip x@( 𝕀 lo hi )
-  | lo == fromInteger 0 && hi == fromInteger 0
+  | lo == 0 && hi == 0
   = [ 𝕀 negInf posInf ]
-  | lo >= fromInteger 0 || hi <= fromInteger 0
+  | lo >= 0 || hi <= 0
   = [ recip x ]
   | otherwise
-  = [ 𝕀 negInf ( recip lo ), 𝕀 ( recip hi ) posInf ]
+  = [ recip $ 𝕀 lo -0, recip $ 𝕀 0 hi ]
   where
-    negInf = fromInteger (-1) / fromInteger 0
-    posInf = fromInteger   1  / fromInteger 0
-{-# SPECIALISE extendedRecip :: 𝕀 Double -> [ 𝕀 Double ] #-}
+    negInf, posInf :: Double
+    negInf = -1 / 0
+    posInf =  1 / 0
 
 -------------------------------------------------------------------------------
 -- Lattices.

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

@@ -107,12 +107,12 @@ instance Prelude.Fractional ( 𝕀 Double ) where
   fromRational r =
     let q = Prelude.fromRational r
     in 𝕀 q q
-  recip (𝕀 lo hi)
+  recip ( 𝕀 lo hi )
 -- #ifdef ASSERTS
     | lo == 0
-    = 𝕀 ( fst $ divI 1 hi ) ( 1 / 0 )
+    = 𝕀 ( fst $ divI 1 hi ) ( 1 Prelude./ 0 )
     | hi == 0
-    = 𝕀 ( -1 / 0 ) ( snd $ divI 1 lo )
+    = 𝕀 ( -1 Prelude./ 0 ) ( snd $ divI 1 lo )
     | lo > 0 || hi < 0
 -- #endif
     = 𝕀 ( fst $ divI 1 hi ) ( snd $ divI 1 lo )
@@ -120,7 +120,7 @@ instance Prelude.Fractional ( 𝕀 Double ) where
     | otherwise
     = error "BAD interval recip; should use extendedRecip instead"
 -- #endif
-  p / q = p * recip q
+  p / q = p * Prelude.recip q
 
 instance Floating ( 𝕀 Double ) where
   sqrt = withHW Prelude.sqrt

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

@@ -126,6 +126,7 @@ instance Num a => Ring ( ViaPrelude a ) where
 
 instance Fractional a => Field ( ViaPrelude a ) where
   fromRational = coerce $ Prelude.fromRational @a
+  recip = coerce $ Prelude.recip @a
   (/) = coerce $ (Prelude./) @a
 
 instance Prelude.Floating a => Floating ( ViaPrelude a ) where