From cebfeb0b7ac04aa440f5f47699cd92c6cee3c90b Mon Sep 17 00:00:00 2001 From: sheaf Date: Wed, 6 Mar 2024 16:07:21 +0100 Subject: [PATCH] Fix some more issues with interval recip --- brush-strokes/src/cusps/bench/Main.hs | 44 ++++++++----------- .../Math/Bezier/Stroke/EnvelopeEquation.hs | 2 + brush-strokes/src/lib/Math/Interval.hs | 17 ++++--- .../src/lib/Math/Interval/Internal.hs | 8 ++-- brush-strokes/src/lib/Math/Ring.hs | 1 + 5 files changed, 34 insertions(+), 38 deletions(-) diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs index 36e793c..fc804e0 100644 --- a/brush-strokes/src/cusps/bench/Main.hs +++ b/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 --} \ No newline at end of file +( b1 - a12 * x2 ) `extendedDivide` a11 +-} diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs b/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs index d5a8adb..fcbb7d6 100644 --- a/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs +++ b/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: diff --git a/brush-strokes/src/lib/Math/Interval.hs b/brush-strokes/src/lib/Math/Interval.hs index 70e6232..df141a6 100644 --- a/brush-strokes/src/lib/Math/Interval.hs +++ b/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. diff --git a/brush-strokes/src/lib/Math/Interval/Internal.hs b/brush-strokes/src/lib/Math/Interval/Internal.hs index d58e8e6..026e090 100644 --- a/brush-strokes/src/lib/Math/Interval/Internal.hs +++ b/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 diff --git a/brush-strokes/src/lib/Math/Ring.hs b/brush-strokes/src/lib/Math/Ring.hs index ae3773a..504713b 100644 --- a/brush-strokes/src/lib/Math/Ring.hs +++ b/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