From 52420a1169ea8bb1de09ace7b7cc535f8e66c0dc Mon Sep 17 00:00:00 2001 From: sheaf Date: Sun, 12 Mar 2023 19:15:58 +0100 Subject: [PATCH] experiment: FMA backend for interval arithmetic also includes the start of a way to observe which equations are being solved, which should help with improving the performance of the interval Newton method --- MetaBrush.cabal | 35 ++- cabal.project | 18 +- src/app/MetaBrush/Application.hs | 4 +- src/app/MetaBrush/Render/Document.hs | 4 +- src/cusps/Main.hs | 209 ++++++++++++++ src/metabrushes/MetaBrush/Asset/Brushes.hs | 13 +- src/metabrushes/MetaBrush/Brush.hs | 23 +- src/metabrushes/MetaBrush/Records.hs | 12 +- src/splines/Math/Algebra/Dual.hs | 2 +- src/splines/Math/Algebra/Dual/Internal.hs | 2 +- src/splines/Math/Bezier/Stroke.hs | 267 +++++++++--------- .../Math/Bezier/Stroke/EnvelopeEquation.hs | 28 +- src/splines/Math/Differentiable.hs | 60 ++-- src/splines/Math/Interval.hs | 60 ++-- src/splines/Math/Interval/FMA.hs | 152 ++++++++++ src/splines/Math/Interval/Internal.hs | 163 +++++++++-- 16 files changed, 770 insertions(+), 282 deletions(-) create mode 100644 src/cusps/Main.hs create mode 100644 src/splines/Math/Interval/FMA.hs diff --git a/MetaBrush.cabal b/MetaBrush.cabal index 487af6a..dbe18e8 100644 --- a/MetaBrush.cabal +++ b/MetaBrush.cabal @@ -29,6 +29,11 @@ flag asserts default: False manual: True +flag use-fma + description: Use fused-muliply add instructions to implement interval arithmetic. + default: False + manual: True + common common build-depends: @@ -53,7 +58,7 @@ common common , rounded-hw ^>= 0.3 , transformers - ^>= 0.5.6.2 + >= 0.5.6.2 && < 0.7 default-extensions: BangPatterns @@ -109,6 +114,14 @@ common common cpp-options: -DASSERTS + if flag(use-fma) + cpp-options: + -DUSE_FMA + ghc-options: + -mfma + if impl(ghc < 9.7) + buildable: False + autogen-modules: Paths_MetaBrush @@ -201,6 +214,10 @@ library splines , Math.Module.Internal , TH.Utils + if flag(use-fma) + other-modules: + Math.Interval.FMA + build-depends: bifunctors >= 5.5.4 && < 5.6 @@ -248,6 +265,19 @@ library metabrushes , bytestring >= 0.10.10.0 && < 0.12 +executable test-cusp-isolation + + import: + common + + hs-source-dirs: + src/cusps + + build-depends: + splines + + main-is: + Main.hs executable convert-metafont @@ -261,8 +291,7 @@ executable convert-metafont Haskell2010 main-is: - - Main.hs + Main.hs other-modules: MetaBrush.MetaFont.Convert diff --git a/cabal.project b/cabal.project index f40c56d..4dcbede 100644 --- a/cabal.project +++ b/cabal.project @@ -11,6 +11,10 @@ source-repository-package location: https://github.com/haskell-waargonaut/waargonaut tag: 5f838582a8c5aae1a198ecd4958729e53a6b03cf +allow-newer: + *:base, *:ghc, *:ghc-prim, *:template-haskell, + *:text + ------------- -- GHC 9.2 -- ------------- @@ -18,11 +22,8 @@ source-repository-package allow-newer: digit:lens, hedgehog:resourcet, - hoist-error:base, natural:lens, - prim-instances:base, - records-sop:ghc-prim, - waargonaut:base, waargonaut:lens, + waargonaut:lens, waargonaut:records-sop, waargonaut:semigroups, waargonaut:text, waargonaut:vector, waargonaut:witherable, @@ -54,9 +55,6 @@ source-repository-package -- GHC 9.6 -- ------------- -allow-newer: - attoparsec:ghc-prim, - generics-sop:base, generics-sop:ghc-prim, - hw-prim:ghc-prim, - integer-logarithms:ghc-prim, - vector-stream:base, vector-stream:ghc-prim, +------------- +-- GHC 9.8 -- +------------- diff --git a/src/app/MetaBrush/Application.hs b/src/app/MetaBrush/Application.hs index ac1b805..3a2697c 100644 --- a/src/app/MetaBrush/Application.hs +++ b/src/app/MetaBrush/Application.hs @@ -170,8 +170,8 @@ runApplication application = do { splineStart = mkPoint ( ℝ2 10 -20 ) 2 1 0 , splineCurves = OpenCurves $ Seq.fromList [ LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 10 10 ) 10 5 ( pi / 4 ) ), curveData = invalidateCache undefined } - , LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10 10 ) 8 5 ( pi / 4 ) ), curveData = invalidateCache undefined } - , LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10 -20 ) 10 7 ( pi / 2 ) ), curveData = invalidateCache undefined } + --, LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10 10 ) 8 5 ( pi / 4 ) ), curveData = invalidateCache undefined } + --, LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10 -20 ) 10 7 ( pi / 2 ) ), curveData = invalidateCache undefined } ] } } diff --git a/src/app/MetaBrush/Render/Document.hs b/src/app/MetaBrush/Render/Document.hs index edea394..d23bead 100644 --- a/src/app/MetaBrush/Render/Document.hs +++ b/src/app/MetaBrush/Render/Document.hs @@ -75,8 +75,6 @@ import Math.Bezier.Stroke ( Cusp(..), CachedStroke(..), invalidateCache , computeStrokeOutline ) -import Math.Interval - ( Extent(Point) ) import Math.Linear ( ℝ(..), T(..) ) import Math.Module @@ -310,7 +308,7 @@ strokeRenderData fitParams StrokeWithOutlineRenderData { strokeDataSpline = spline , strokeOutlineData = ( outline, fitPts, cusps ) - , strokeBrushFunction = fun @Double ( brushFn @Point proxy# id ) + , strokeBrushFunction = fun @Double ( brushFn @2 @() proxy# id ) . embedUsedParams . toUsedParams } diff --git a/src/cusps/Main.hs b/src/cusps/Main.hs new file mode 100644 index 0000000..fb72764 --- /dev/null +++ b/src/cusps/Main.hs @@ -0,0 +1,209 @@ +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} + +module Main where + +-- base +import Prelude hiding ( Num(..), Fractional(..), Floating(..), (^) ) +import qualified Prelude ( Num(..), Fractional(..) ) +import Data.Coerce + ( Coercible ) +import GHC.Generics + ( Generic ) + +-- acts +import Data.Act + ( Torsor ) + +-- containers +import Data.Map + ( Map ) +import qualified Data.Map as Map + ( (!) ) +import qualified Data.Sequence as Seq + ( index ) + +-- generic-lens +import Data.Generics.Product.Typed + ( HasType ) + +-- MetaBrush +import Math.Algebra.Dual +import Math.Bezier.Spline +import Math.Bezier.Stroke + ( brushStrokeData, pathAndUsedParams ) +import Math.Bezier.Stroke.EnvelopeEquation + ( StrokeDatum(..) ) +import Math.Differentiable + ( I ) +import Math.Interval +import Math.Linear + ( ℝ(..), T(..) ) +import Math.Module + ( Module, Cross ) +import Math.Ring + ( AbelianGroup(..), Ring(..), Field(..), Transcendental(..) + , ViaPrelude(..) + ) + +-------------------------------------------------------------------------------- + +data PointData params + = PointData + { pointCoords :: !( ℝ 2 ) + , brushParams :: !params + } + deriving stock ( Show, Generic ) + + +{- + +newtype C k u v = D { runD :: u -> D k u v } + +type instance D k ( ℝ n ) = Dk𝔸n +type instance D k ( 𝕀 v ) = D k v + +e.g. + +data D2𝔸2 v = + D22 { _D22_v :: !v + , _D22_dx, _D22_dy :: !( T v ) + , _D22_dxdx, _D22_dxdy, _D22_dydy :: !( T v ) + } + +-} + + +outlineFunction + :: forall i brushParams + . ( Show brushParams + , D 1 ( I i ( ℝ 2 ) ) ~ D 1 ( ℝ 2 ) + , D 2 ( I i ( ℝ 2 ) ) ~ D 2 ( ℝ 2 ) + , D 3 ( I i ( ℝ 1 ) ) ~ D 3 ( ℝ 1 ) + , D 3 ( I i ( ℝ 2 ) ) ~ D 3 ( ℝ 2 ) + , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) + , HasType ( ℝ 2 ) ( PointData brushParams ) + , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) + , Module ( I i Double ) ( T ( I i brushParams ) ) + , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) ) + , Torsor ( T ( I i brushParams ) ) ( I i brushParams ) + , HasChainRule ( I i Double ) 3 ( I i ( ℝ 1 ) ) + , HasChainRule ( I i Double ) 3 ( I i brushParams ) + , Traversable ( D 3 brushParams ) + , Traversable ( D 3 ( I i brushParams ) ) + ) + => ( forall a. a -> I i a ) + -> C 3 ( I i brushParams ) ( Spline Closed () ( I i ( ℝ 2 ) ) ) + -- ^ brush shape + -> Int -- ^ brush segment index to consider + -> PointData brushParams -> Curve Open () ( PointData brushParams ) + -- ^ stroke path + -> ( I i ( ℝ 1 ) -> I i ( ℝ 1 ) -> StrokeDatum 3 i ) +outlineFunction single brush brush_index sp0 crv = strokeData + + where + path :: C 3 ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) + params :: C 3 ( I i ( ℝ 1 ) ) ( I i brushParams ) + (path, params) = + pathAndUsedParams @3 @i single brushParams sp0 crv + + strokeData :: I i ( ℝ 1 ) -> I i ( ℝ 1 ) -> StrokeDatum 3 i + strokeData = + fmap ( `Seq.index` brush_index ) $ + brushStrokeData @3 @brushParams @i path params brush + +main :: IO () +main = return () + +-------------------------------------------------------------------------------- +-- EDSL for inspection + +type instance I AA a = AA a + +data AAVal a where + Val :: a -> AAVal a + Var :: String -> AAVal Double + +instance Show a => Show ( AAVal a ) where + show (Val v) = show v + show (Var v) = show v + +data AA a where + Pt :: AAVal a -> AA a + Ival :: AAVal a -> AAVal a -> AA a + (:+:) :: AA Double -> AA Double -> AA Double + (:-:) :: AA Double -> AA Double -> AA Double + Negate :: AA Double -> AA Double + (:*:) :: AA Double -> AA Double -> AA Double + (:^:) :: AA Double -> Word -> AA Double + Recip :: AA Double -> AA Double + Cos :: AA Double -> AA Double + Sin :: AA Double -> AA Double + Pi :: AA Double + +instance Show a => Show (AA a) where + showsPrec prec = \case + Pt v -> showString "[" . showsPrec 0 v . showString "]" + Ival lo hi -> showString "[" . showsPrec 0 lo . showString "," . showsPrec 0 hi . showString "]" + iv1 :+: iv2 -> showParen (prec > 6) $ showsPrec 6 iv1 . showString " + " . showsPrec 7 iv2 + iv1 :-: iv2 -> showParen (prec > 6) $ showsPrec 6 iv1 . showString " - " . showsPrec 7 iv2 + iv1 :*: iv2 -> showParen (prec > 7) $ showsPrec 7 iv1 . showString " - " . showsPrec 8 iv2 + iv :^: n -> showParen (prec > 8) $ showsPrec 9 iv . showString " ^ " . showsPrec 8 n + Negate iv -> showParen (prec > 10) $ showString "negate " . showsPrec 11 iv + Recip iv -> showParen (prec > 10) $ showString "recip " . showsPrec 11 iv + Cos iv -> showParen (prec > 10) $ showString "cos " . showsPrec 11 iv + Sin iv -> showParen (prec > 10) $ showString "sin " . showsPrec 11 iv + Pi -> showString "pi" + +infixl 6 :+: +infixl 6 :-: +infixl 7 :*: +infixr 8 :^: + +val :: Map String Double -> AAVal a -> a +val _ ( Val a ) = a +val vars ( Var v ) = vars Map.! v + +eval :: Map String Double -> AA a -> 𝕀 a +eval vars = \ case + Pt v -> let f = val vars v in 𝕀 f f + Ival lo hi -> 𝕀 (val vars lo) (val vars hi) + iv1 :+: iv2 -> eval vars iv1 + eval vars iv2 + iv1 :-: iv2 -> eval vars iv1 - eval vars iv2 + Negate iv -> negate $ eval vars iv + iv1 :*: iv2 -> eval vars iv1 * eval vars iv2 + iv :^: n -> eval vars iv ^ n + Recip iv -> recip $ eval vars iv + Cos iv -> cos $ eval vars iv + Sin iv -> sin $ eval vars iv + Pi -> pi + +instance Prelude.Num ( AA Double ) where + (+) = (:+:) + (-) = (:-:) + negate = Negate + (*) = (:*:) + abs = error "No abs for abstract intervals" + signum = error "No signum for abstract intervals" + fromInteger = Pt . Val . fromInteger + +instance Prelude.Fractional ( AA Double ) where + recip = Recip + (/) = error "No division for abstract intervals" + fromRational = Pt . Val . fromRational + +instance Ring ( AA Double ) where + (*) = (:*:) + (^) = (:^:) + +deriving via ViaPrelude ( AA Double ) + instance AbelianGroup ( AA Double ) +deriving via ViaPrelude ( AA Double ) + instance AbelianGroup ( T ( AA Double ) ) +deriving via ViaPrelude ( AA Double ) + instance Field ( AA Double ) + +instance Transcendental ( AA Double ) where + pi = Pi + cos = Cos + sin = Sin diff --git a/src/metabrushes/MetaBrush/Asset/Brushes.hs b/src/metabrushes/MetaBrush/Asset/Brushes.hs index 31bc2fa..b10e50f 100644 --- a/src/metabrushes/MetaBrush/Asset/Brushes.hs +++ b/src/metabrushes/MetaBrush/Asset/Brushes.hs @@ -1,5 +1,6 @@ {-# LANGUAGE AllowAmbiguousTypes #-} {-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE PolyKinds #-} {-# LANGUAGE ScopedTypeVariables #-} module MetaBrush.Asset.Brushes @@ -32,8 +33,8 @@ import qualified Data.HashMap.Strict as HashMap -- MetaBrush import Math.Algebra.Dual import Math.Bezier.Spline -import Math.Interval - +import Math.Differentiable + ( I ) import Math.Linear import Math.Module ( Module((^+^), (*^)) ) @@ -81,7 +82,7 @@ ellipse = BrushData "ellipse" ( WithParams deflts ellipseBrush ) -------------------------------------------------------------------------------- -- Differentiable brushes. -circleSpline :: forall i k u v +circleSpline :: forall {t} (i :: t) k u v . Applicative ( D k ( I i u ) ) => ( Double -> Double -> D k ( I i u ) ( I i v ) ) -> D k ( I i u ) ( Spline 'Closed () ( I i v ) ) @@ -98,7 +99,7 @@ circleSpline p = sequenceA $ Bezier3To ( p κ -1 ) ( p 1 -κ ) BackToStart () {-# INLINE circleSpline #-} -circleBrush :: forall i k irec +circleBrush :: forall {t} (i :: t) k irec . ( irec ~ I i ( Record CircleBrushFields ) , Module ( D k irec ( I i Double ) ) @@ -113,7 +114,7 @@ circleBrush :: forall i k irec -> C k irec ( Spline 'Closed () ( I i ( ℝ 2 ) ) ) circleBrush _ mkI = D \ params -> - let r :: D k irec( I i Double ) + let r :: D k irec ( I i Double ) r = runD ( var @_ @k ( Fin 1 ) ) params mkPt :: Double -> Double -> D k irec ( I i ( ℝ 2 ) ) mkPt ( kon -> x ) ( kon -> y ) @@ -128,7 +129,7 @@ circleBrush _ mkI = kon = konst @( I i Double ) @k @irec . mkI {-# INLINEABLE circleBrush #-} -ellipseBrush :: forall i k irec +ellipseBrush :: forall {t} (i :: t) k irec . ( irec ~ I i ( Record EllipseBrushFields ) , Module ( D k irec ( I i Double ) ) diff --git a/src/metabrushes/MetaBrush/Brush.hs b/src/metabrushes/MetaBrush/Brush.hs index 8ea2d0c..d9618e9 100644 --- a/src/metabrushes/MetaBrush/Brush.hs +++ b/src/metabrushes/MetaBrush/Brush.hs @@ -1,4 +1,5 @@ {-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE PolyKinds #-} {-# LANGUAGE QuantifiedConstraints #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE UndecidableInstances #-} @@ -44,9 +45,9 @@ import Math.Algebra.Dual import Math.Bezier.Spline ( SplineType(Closed), Spline ) import Math.Differentiable - ( DiffInterp, ExtentOrder ) + ( DiffInterp, I ) import Math.Interval - ( type I, Extent(Point, Interval) ) + ( 𝕀 ) import Math.Linear import MetaBrush.Records ( KnownSymbols, Length, Record ) @@ -62,11 +63,11 @@ data WithParams params f = WithParams { defaultParams :: params , withParams - :: forall i - . ( DiffInterp i params ) + :: forall {t} k (i :: t) + . ( DiffInterp k i params ) => Proxy# i -> ( forall a. a -> I i a ) - -> C ( ExtentOrder i ) ( I i params ) ( f ( I i ( ℝ 2 ) ) ) + -> C k ( I i params ) ( f ( I i ( ℝ 2 ) ) ) } -------------------------------------------------------------------------------- @@ -82,8 +83,8 @@ data Brush brushFields where :: forall brushFields . ( KnownSymbols brushFields, Typeable brushFields , Representable Double ( ℝ ( Length brushFields ) ) - , DiffInterp Point ( ℝ ( Length brushFields ) ) - , DiffInterp Interval ( ℝ ( Length brushFields ) ) + , DiffInterp 2 () ( ℝ ( Length brushFields ) ) + , DiffInterp 3 𝕀 ( ℝ ( Length brushFields ) ) ) => { brushName :: !Text , brushFunction :: BrushFunction brushFields @@ -118,8 +119,8 @@ class ( KnownSymbols pointFields, Typeable pointFields , Show ( Record pointFields ) , NFData ( Record pointFields ) , Representable Double ( ℝ ( Length pointFields ) ) - , DiffInterp Point ( ℝ ( Length pointFields ) ) - , DiffInterp Interval ( ℝ ( Length pointFields ) ) + , DiffInterp 2 () ( ℝ ( Length pointFields ) ) + , DiffInterp 3 𝕀 ( ℝ ( Length pointFields ) ) ) => PointFields pointFields where { } instance ( KnownSymbols pointFields, Typeable pointFields @@ -127,8 +128,8 @@ instance ( KnownSymbols pointFields, Typeable pointFields , Show ( Record pointFields ) , NFData ( Record pointFields ) , Representable Double ( ℝ ( Length pointFields ) ) - , DiffInterp Point ( ℝ ( Length pointFields ) ) - , DiffInterp Interval ( ℝ ( Length pointFields ) ) + , DiffInterp 2 () ( ℝ ( Length pointFields ) ) + , DiffInterp 3 𝕀 ( ℝ ( Length pointFields ) ) ) => PointFields pointFields where { } diff --git a/src/metabrushes/MetaBrush/Records.hs b/src/metabrushes/MetaBrush/Records.hs index 49711fc..6017317 100644 --- a/src/metabrushes/MetaBrush/Records.hs +++ b/src/metabrushes/MetaBrush/Records.hs @@ -181,8 +181,8 @@ intersect :: forall r1 r2 l1 l2 , l1 ~ Length r1, l2 ~ Length r2 , Representable Double ( ℝ l1 ) , Representable Double ( ℝ l2 ) - , Differentiable 'Point ( ℝ l2 ) - , Differentiable 'Interval ( ℝ l2 ) + , Differentiable 2 () ( ℝ l2 ) + , Differentiable 3 𝕀 ( ℝ l2 ) ) => Intersection r1 r2 intersect @@ -206,8 +206,8 @@ data Intersection r1 r2 where . ( l12 ~ Length r1r2 , KnownSymbols r1r2 , Representable Double ( ℝ l12 ) - , Differentiable 'Point ( ℝ l12 ) - , Differentiable 'Interval ( ℝ l12 ) + , Differentiable 2 () ( ℝ l12 ) + , Differentiable 3 𝕀 ( ℝ l12 ) ) => { project :: Record r1 -> Record r1r2 -- ^ project out fields present in both rows @@ -226,8 +226,8 @@ doIntersection => ( forall r1r2 l12. ( r1r2 ~ Intersect r1 r2 , KnownSymbols r1r2, l12 ~ Length r1r2 - , Differentiable 'Point ( ℝ l12 ) - , Differentiable 'Interval ( ℝ l12 ) + , Differentiable 2 () ( ℝ l12 ) + , Differentiable 3 𝕀 ( ℝ l12 ) , Representable Double ( ℝ l12 ) ) => Proxy# r1r2 -> Vec l12 ( Fin l1 ) -> Vec l12 ( Fin l2 ) -> kont ) diff --git a/src/splines/Math/Algebra/Dual.hs b/src/splines/Math/Algebra/Dual.hs index 3f32195..64a2402 100644 --- a/src/splines/Math/Algebra/Dual.hs +++ b/src/splines/Math/Algebra/Dual.hs @@ -6,7 +6,7 @@ {-# OPTIONS_GHC -Wno-orphans -O2 #-} -{-# OPTIONS_GHC -ddump-splices -ddump-simpl -ddump-to-file -dno-typeable-binds +{-# OPTIONS_GHC -ddump-splices -ddump-simpl -ddump-cmm -ddump-to-file -dno-typeable-binds -dsuppress-unfoldings -dsuppress-coercions #-} module Math.Algebra.Dual diff --git a/src/splines/Math/Algebra/Dual/Internal.hs b/src/splines/Math/Algebra/Dual/Internal.hs index 723ab14..32c4814 100644 --- a/src/splines/Math/Algebra/Dual/Internal.hs +++ b/src/splines/Math/Algebra/Dual/Internal.hs @@ -178,7 +178,7 @@ data D2𝔸4 v = deriving Applicative via Generically1 D2𝔸4 --- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^3 \). +-- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^4 \). data D3𝔸4 v = D34 { _D34_v :: !v , _D34_dx, _D34_dy, _D34_dz, _D34_dw :: !( T v ) diff --git a/src/splines/Math/Bezier/Stroke.hs b/src/splines/Math/Bezier/Stroke.hs index 93da02a..9f20ba0 100644 --- a/src/splines/Math/Bezier/Stroke.hs +++ b/src/splines/Math/Bezier/Stroke.hs @@ -1,4 +1,5 @@ {-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE PolyKinds #-} {-# LANGUAGE QuantifiedConstraints #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE UndecidableInstances #-} @@ -13,6 +14,7 @@ module Math.Bezier.Stroke , brushStroke, envelopeEquation , line, bezier2, bezier3 + , brushStrokeData, pathAndUsedParams ) where @@ -45,9 +47,8 @@ import Data.Maybe import Data.Semigroup ( sconcat ) import GHC.Exts - ( newMutVar#, runRW# + ( newMutVar#, runRW#, inline , Proxy#, proxy# - , inline ) import GHC.STRef ( STRef(..), readSTRef, writeSTRef ) @@ -88,13 +89,6 @@ import Data.Generics.Internal.VL import qualified Control.Parallel.Strategies as Strats ( rdeepseq, parTuple2, using ) --- rounded-hw -import Numeric.Rounded.Hardware - ( Rounded(..) ) -import Numeric.Rounded.Hardware.Interval.NonEmpty - ( Interval(..) ) -import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval - -- transformers import Control.Monad.Trans.Class ( lift ) @@ -123,9 +117,7 @@ import Math.Bezier.Spline import qualified Math.Bezier.Quadratic as Quadratic import Math.Bezier.Stroke.EnvelopeEquation import Math.Differentiable - ( Differentiable, DiffInterp - , type ExtentOrder - ) + ( Differentiable, DiffInterp, I ) import Math.Epsilon ( epsilon, nearZero ) import Math.Interval @@ -236,13 +228,13 @@ computeStrokeOutline :: -- Differentiability. , Interpolatable Double usedParams , Interpolatable ( 𝕀 Double ) ( 𝕀 usedParams ) - , DiffInterp 'Point brushParams - , DiffInterp 'Interval brushParams - , HasChainRule Double ( ExtentOrder 'Point ) usedParams - , HasChainRule ( 𝕀 Double ) ( ExtentOrder 'Interval ) ( 𝕀 usedParams ) - , HasChainRule Double ( ExtentOrder 'Point ) brushParams - , HasChainRule ( 𝕀 Double ) ( ExtentOrder 'Interval ) ( 𝕀 brushParams ) - , Traversable ( D ( ExtentOrder 'Point ) brushParams ) + , DiffInterp 2 () brushParams + , DiffInterp 3 𝕀 brushParams + , HasChainRule Double 2 usedParams + , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 usedParams ) + , HasChainRule Double 2 brushParams + , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 brushParams ) + , Traversable ( D 2 brushParams ) -- Debugging. , Show ptData, Show brushParams @@ -251,13 +243,12 @@ computeStrokeOutline :: => FitParameters -> ( ptData -> usedParams ) -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing - -> ( forall i - . DiffInterp i brushParams + -> ( forall {t} k (i :: t) + . DiffInterp k i brushParams => Proxy# i -> ( forall a. a -> I i a ) - -> C ( ExtentOrder i ) - ( I i brushParams ) - ( Spline Closed () ( I i ( ℝ 2 ) ) ) + -> C k ( I i brushParams ) + ( Spline Closed () ( I i ( ℝ 2 ) ) ) ) -> Spline clo crvData ptData -> ST s @@ -381,7 +372,7 @@ computeStrokeOutline fitParams ptParams toBrushParams brushFn spline@( Spline { outlineInfo p0 crv :<| go ( openCurveEnd crv ) crvs brushShape :: ptData -> SplinePts Closed - brushShape pt = fun @Double ( brushFn @Point proxy# id ) $ toBrushParams $ ptParams pt + brushShape pt = fun @Double ( brushFn @2 @() proxy# id ) $ toBrushParams $ ptParams pt updateSpline :: ( T ( ℝ 2 ), T ( ℝ 2 ), T ( ℝ 2 ) ) -> ST s OutlineData updateSpline ( lastTgt, lastTgtFwd, lastTgtBwd ) @@ -514,87 +505,61 @@ outlineFunction -- Differentiability. , Interpolatable Double usedParams , Interpolatable ( 𝕀 Double ) ( 𝕀 usedParams ) - , DiffInterp 'Point brushParams - , DiffInterp 'Interval brushParams - , HasChainRule Double ( ExtentOrder 'Point ) usedParams - , HasChainRule ( 𝕀 Double ) ( ExtentOrder 'Interval ) ( 𝕀 usedParams ) - , HasChainRule Double ( ExtentOrder 'Point ) brushParams - , HasChainRule ( 𝕀 Double ) ( ExtentOrder 'Interval ) ( 𝕀 brushParams ) - , Traversable ( D ( ExtentOrder 'Point ) brushParams ) + , DiffInterp 2 () brushParams + , DiffInterp 3 𝕀 brushParams + , HasChainRule Double 2 usedParams + , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 usedParams ) + , HasChainRule Double 2 brushParams + , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 brushParams ) + , Traversable ( D 2 brushParams ) -- Debugging. , Show ptData, Show brushParams ) => ( ptData -> usedParams ) -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing - -> ( forall i - . DiffInterp i brushParams + -> ( forall {t} k (i :: t) + . DiffInterp k i brushParams => Proxy# i -> ( forall a. a -> I i a ) - -> C ( ExtentOrder i ) - ( I i brushParams ) - ( Spline Closed () ( I i ( ℝ 2 ) ) ) + -> C k ( I i brushParams ) + ( Spline Closed () ( I i ( ℝ 2 ) ) ) ) -> ptData -> Curve Open crvData ptData -> OutlineInfo outlineFunction ptParams toBrushParams brushFromParams = \ sp0 crv -> let - pathAndUsedParams :: forall i k arr - . ( k ~ ExtentOrder i, CurveOrder k - , arr ~ C k - , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) - , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) - , Module ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) - , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) ) - , Module ( I i Double ) ( T ( I i usedParams ) ) - , Torsor ( T ( I i usedParams ) ) ( I i usedParams ) - ) - => ( forall a. a -> I i a ) - -> ( I i ( ℝ 1 ) `arr` I i ( ℝ 2 ), I i ( ℝ 1 ) `arr` I i usedParams ) - pathAndUsedParams toI = - case crv of - LineTo { curveEnd = NextPoint sp1 } - | let seg = Segment sp0 sp1 - -> ( line @k @i ( fmap ( toI . coords ) seg ) - , line @k @i ( fmap ( toI . ptParams ) seg ) ) - Bezier2To { controlPoint = sp1, curveEnd = NextPoint sp2 } - | let bez2 = Quadratic.Bezier sp0 sp1 sp2 - -> ( bezier2 @k @i ( fmap ( toI . coords ) bez2 ) - , bezier2 @k @i ( fmap ( toI . ptParams ) bez2 ) ) - Bezier3To { controlPoint1 = sp1, controlPoint2 = sp2, curveEnd = NextPoint sp3 } - | let bez3 = Cubic.Bezier sp0 sp1 sp2 sp3 - -> ( bezier3 @k @i ( fmap ( toI . coords ) bez3 ) - , bezier3 @k @i ( fmap ( toI . ptParams ) bez3 ) ) - usedParams :: C ( ExtentOrder 'Point ) ( ℝ 1 ) usedParams - path :: C ( ExtentOrder 'Point ) ( ℝ 1 ) ( ℝ 2 ) - ( path, usedParams ) = pathAndUsedParams @Point id + usedParams :: C 2 ( ℝ 1 ) usedParams + path :: C 2 ( ℝ 1 ) ( ℝ 2 ) + ( path, usedParams ) + = pathAndUsedParams @2 @() id ptParams sp0 crv curves :: ℝ 1 -- t - -> Seq ( ℝ 1 {- s -} -> StrokeDatum Point ) + -> Seq ( ℝ 1 {- s -} -> StrokeDatum 2 () ) curves = - brushStrokeData @Point @( ExtentOrder 'Point ) @brushParams + brushStrokeData @2 @brushParams path - ( chainRule @Double @( ExtentOrder 'Point ) + ( chainRule @Double @2 usedParams ( linear toBrushParams ) ) - ( brushFromParams @Point proxy# id ) + ( brushFromParams @2 @() proxy# id ) curvesI :: 𝕀ℝ 1 -- t - -> Seq ( 𝕀ℝ 1 {- s -} -> StrokeDatum 'Interval ) - curvesI = brushStrokeData @'Interval @( ExtentOrder 'Interval ) @brushParams + -> Seq ( 𝕀ℝ 1 {- s -} -> StrokeDatum 3 𝕀 ) + curvesI = brushStrokeData @3 @brushParams pathI - ( chainRule @( 𝕀 Double ) @( ExtentOrder 'Interval ) + ( chainRule @( 𝕀 Double ) @3 usedParamsI ( linear ( nonDecreasing toBrushParams ) ) ) - ( brushFromParams @'Interval proxy# singleton ) + ( brushFromParams @3 @𝕀 proxy# singleton ) - usedParamsI :: C ( ExtentOrder 'Interval ) ( 𝕀ℝ 1 ) ( 𝕀 usedParams ) - pathI :: C ( ExtentOrder 'Interval ) ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 ) - ( pathI, usedParamsI ) = pathAndUsedParams @'Interval singleton + usedParamsI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀 usedParams ) + pathI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 ) + ( pathI, usedParamsI ) = pathAndUsedParams @3 @𝕀 singleton ptParams sp0 crv fwdBwd :: OutlineFn fwdBwd t @@ -608,10 +573,10 @@ outlineFunction ptParams toBrushParams brushFromParams = \ sp0 crv -> D21 path_t path'_t _ = runD path t D21 params_t _ _ = runD usedParams t brush_t = value @Double @2 @brushParams - $ runD ( brushFromParams @Point proxy# id ) + $ runD ( brushFromParams @2 @() proxy# id ) $ toBrushParams params_t - ( newtDunno, newtSols ) = intervalNewtonGS InverseMidJacobian 0.0001 curvesI + ( newtDunno, newtSols ) = intervalNewtonGS InverseMidJacobian 1e-7 curvesI in --trace -- ( unlines $ @@ -631,6 +596,37 @@ outlineFunction ptParams toBrushParams brushFromParams = \ sp0 crv -> } {-# INLINEABLE outlineFunction #-} +pathAndUsedParams :: forall k i arr crvData ptData usedParams + . ( HasType ( ℝ 2 ) ptData + , CurveOrder k + , arr ~ C k + , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) + , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) + , Module ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) + , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) ) + , Module ( I i Double ) ( T ( I i usedParams ) ) + , Torsor ( T ( I i usedParams ) ) ( I i usedParams ) + ) + => ( forall a. a -> I i a ) + -> ( ptData -> usedParams ) + -> ptData + -> Curve Open crvData ptData + -> ( I i ( ℝ 1 ) `arr` I i ( ℝ 2 ), I i ( ℝ 1 ) `arr` I i usedParams ) +pathAndUsedParams toI ptParams sp0 crv = + case crv of + LineTo { curveEnd = NextPoint sp1 } + | let seg = Segment sp0 sp1 + -> ( line @k @i ( fmap ( toI . coords ) seg ) + , line @k @i ( fmap ( toI . ptParams ) seg ) ) + Bezier2To { controlPoint = sp1, curveEnd = NextPoint sp2 } + | let bez2 = Quadratic.Bezier sp0 sp1 sp2 + -> ( bezier2 @k @i ( fmap ( toI . coords ) bez2 ) + , bezier2 @k @i ( fmap ( toI . ptParams ) bez2 ) ) + Bezier3To { controlPoint1 = sp1, controlPoint2 = sp2, curveEnd = NextPoint sp3 } + | let bez3 = Cubic.Bezier sp0 sp1 sp2 sp3 + -> ( bezier3 @k @i ( fmap ( toI . coords ) bez3 ) + , bezier3 @k @i ( fmap ( toI . ptParams ) bez3 ) ) + ----------------------------------- -- Various utility functions -- used in the "stroke" function. @@ -924,10 +920,8 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) , offset = T $ Cubic.bezier @( T ( ℝ 2 ) ) bez t } - - -splineCurveFns :: forall i k - . ( k ~ ExtentOrder i, CurveOrder k +splineCurveFns :: forall k i + . ( CurveOrder k , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) , Module ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) ) @@ -958,7 +952,7 @@ solveEnvelopeEquations :: ℝ 1 -- ^ @t@ (for debugging only) -> ℝ 2 -> T ( ℝ 2 ) -> ( Offset, Offset ) - -> Seq ( ℝ 1 -> StrokeDatum Point ) + -> Seq ( ℝ 1 -> StrokeDatum 2 () ) -> ( ( ℝ 2, T ( ℝ 2 ) ), ( ℝ 2, T ( ℝ 2 ) ) ) solveEnvelopeEquations _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData = ( fwdSol, ( bwdPt, -1 *^ bwdTgt ) ) @@ -1020,7 +1014,7 @@ solveEnvelopeEquations _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData | otherwise = i + 1 - sol :: Double -> ( ℝ 1 -> StrokeDatum Point ) -> ( Bool, ℝ 1, ℝ 2, T ( ℝ 2 ) ) + sol :: Double -> ( ℝ 1 -> StrokeDatum 2 () ) -> ( Bool, ℝ 1, ℝ 2, T ( ℝ 2 ) ) sol initialGuess f = let (good, s) = case newtonRaphson maxIters precision domain ( eqn f ) initialGuess of Nothing -> ( False, initialGuess ) @@ -1052,13 +1046,13 @@ solveEnvelopeEquations _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData -- , " ∂E/∂s = " ++ show ee_s -- , " dc/dt = " ++ show dcdt -- ] ) - ( good, ℝ1 s, value @Double @2 @( ℝ 2 ) dstroke, dcdt ) + ( good, ℝ1 s, value @Double @2 @( ℝ 2 ) dstroke, dcdt ) - eqn :: ( ℝ 1 -> StrokeDatum Point ) -> ( Double -> ( Double, Double ) ) + eqn :: ( ℝ 1 -> StrokeDatum 2 () ) -> ( Double -> ( Double, Double ) ) eqn f s = case f ( ℝ1 s ) of - StrokeDatum { ee = D12 ( ℝ1 ee ) _ ( T ( ℝ1 ee_s ) ) } -> - ( ee, ee_s ) + StrokeDatum { ee = D12 ee _ ee_s } -> + coerce ( ee, ee_s ) maxIters :: Word maxIters = 5 --30 @@ -1074,10 +1068,9 @@ instance Applicative ZipSeq where liftA2 f ( ZipSeq xs ) ( ZipSeq ys ) = ZipSeq ( Seq.zipWith f xs ys ) {-# INLINE liftA2 #-} -brushStrokeData :: forall i k brushParams arr - . ( k ~ ExtentOrder i, CurveOrder k, arr ~ C k - , Differentiable i brushParams - , Fractional ( I i Double ) +brushStrokeData :: forall k brushParams i arr + . ( CurveOrder k, arr ~ C k + , Differentiable k i brushParams , HasChainRule ( I i Double ) k ( I i ( ℝ 1 ) ) , Applicative ( D k ( ℝ 1 ) ) @@ -1097,7 +1090,7 @@ brushStrokeData :: forall i k brushParams arr -- ^ brush parameters -> ( I i brushParams `arr` Spline Closed () ( I i ( ℝ 2 ) ) ) -- ^ brush from parameters - -> ( I i ( ℝ 1 ) -> Seq ( I i ( ℝ 1 ) -> StrokeDatum i ) ) + -> ( I i ( ℝ 1 ) -> Seq ( I i ( ℝ 1 ) -> StrokeDatum k i ) ) brushStrokeData path params brush = \ t -> let @@ -1108,7 +1101,7 @@ brushStrokeData path params brush = dbrush_params :: D k ( I i brushParams ) ( Spline Closed () ( I i ( ℝ 2 ) ) ) !dbrush_params = runD brush $ value @( I i Double ) @k @( I i ( ℝ 1 ) ) dparams_t splines :: Seq ( D k ( I i brushParams ) ( I i ( ℝ 1 ) `arr` I i ( ℝ 2 ) ) ) - !splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @i @k ) dbrush_params + !splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @k @i ) dbrush_params dbrushes_t :: Seq ( I i ( ℝ 1 ) -> D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) ) !dbrushes_t = force $ fmap ( uncurryD @k . ( chain @(I i Double) @k dparams_t ) ) splines -- This is the crucial use of the chain rule. @@ -1118,7 +1111,7 @@ brushStrokeData path params brush = mkStrokeDatum :: D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) -> ( I i ( ℝ 1 ) -> D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) ) - -> ( I i ( ℝ 1 ) -> StrokeDatum i ) + -> ( I i ( ℝ 1 ) -> StrokeDatum k i ) mkStrokeDatum dpath_t dbrush_t s = let dbrush_t_s = dbrush_t s dstroke = brushStroke @k dpath_t dbrush_t_s @@ -1200,7 +1193,7 @@ extendedRecip x@( 𝕀 lo hi ) -- | Computes the brush stroke coordinates of a cusp from -- the @(t,s)@ parameter values. -cuspCoords :: ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 'Point ) ) +cuspCoords :: ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) ) -> ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) -> Cusp cuspCoords eqs ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), i, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) ) @@ -1210,8 +1203,8 @@ cuspCoords eqs ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), i, 𝕀 ( ℝ1 s_lo ) ( ℝ1 <- ( eqs ( ℝ1 t_mid ) `Seq.index` i ) ( ℝ1 s_mid ) = Cusp { cuspParameters = ℝ2 t_mid s_mid - , cuspPathCoords = dpath - , cuspStrokeCoords = stroke + , cuspPathCoords = coerce dpath + , cuspStrokeCoords = coerce stroke } where t_mid = 0.5 * ( t_lo + t_hi ) @@ -1228,33 +1221,41 @@ data Preconditioner -- -- Returns @(dunno, sols)@ where @sols@ are boxes that contain a unique solution -- (and to which Newton's method will converge starting from anywhere inside --- the box), and @dunno@ which are small boxes which might or might not +-- the box), and @dunno@ are small boxes which might or might not -- contain solutions. intervalNewtonGS :: Preconditioner -> Double - -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 'Interval ) ) + -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) -> ( [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ], [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] ) intervalNewtonGS precondMethod minWidth eqs = - go [ ( 𝕀 ( ℝ1 0 ) ( ℝ1 1 ), i, 𝕀 ( ℝ1 0 ) ( ℝ1 1 ) ) - | i <- [ 0 .. length ( eqs ( 𝕀 ( ℝ1 0 ) ( ℝ1 1 ) ) ) - 1 ] - ] - [] - [] + go (0,0) + [ ( 𝕀 ( ℝ1 0 ) ( ℝ1 1 ), i, 𝕀 ( ℝ1 0 ) ( ℝ1 1 ) ) + | i <- [ 0 .. length ( eqs ( 𝕀 ( ℝ1 0 ) ( ℝ1 1 ) ) ) - 1 ] + ] + [] + [] where - go :: [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] -- boxes to work on + go :: ( Int, Int ) -- step counts (for debugging) + -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] -- boxes to work on -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] -- too small: don't shrink further -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] -- found solutions -> ( [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ], [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] ) - go [] giveUp sols = ( giveUp, sols ) - go ( cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) ) - , i - , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) ) - ) : cands ) giveUp sols + go ( bis, newt ) [] giveUp sols = + trace ( unlines [ "intervalNewtonGS done" + , " #bisections: " ++ show bis + , " #newtonSteps: " ++ show newt + , " #sols: " ++ show ( length sols ) + , " #dunno: " ++ show ( length giveUp ) ] ) + ( giveUp, sols ) + go ( bis, newt ) ( cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) ) + , i + , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) ) + ) : cands ) giveUp sols -- Box is small: stop processing it. | t_hi - t_lo < minWidth && s_hi - s_lo < minWidth - = go cands ( cand : giveUp ) sols + = go ( bis, newt ) cands ( cand : giveUp ) sols | StrokeDatum { ee = D22 ee _ _ _ _ _ , 𝛿E𝛿sdcdt = D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) } @@ -1262,14 +1263,14 @@ intervalNewtonGS precondMethod minWidth eqs = , StrokeDatum { 𝛿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 - = if | Interval.inf ( ival ee ) < Rounded ( ℝ1 0 ) - , Interval.sup ( ival ee ) > Rounded ( ℝ1 0 ) - , cmpℝ2 (<) ( getRounded ( Interval.inf $ ival f ) ) ( ℝ2 0 0 ) - , cmpℝ2 (>) ( getRounded ( Interval.sup $ ival f ) ) ( ℝ2 0 0 ) + = if | inf ee < ℝ1 0 + , sup ee > ℝ1 0 + , cmpℝ2 (<) ( inf f ) ( ℝ2 0 0 ) + , cmpℝ2 (>) ( sup f ) ( ℝ2 0 0 ) -> let -- Interval Newton method: take one Gauss–Seidel step -- for the equation f'(X) v = - f(x_mid). !( a, b ) = precondition precondMethod - ( f_t_mid, f_s_mid ) + ( f_t, f_s ) ( f_t, f_s ) ( neg f_mid ) !gsGuesses = gaussSeidel a b @@ -1284,7 +1285,7 @@ intervalNewtonGS precondMethod minWidth eqs = -- Newton's method is guaranteed to converge to the unique solution. let !(done, todo) = bimap ( map ( mkGuess . fst ) ) ( map ( mkGuess . fst ) ) $ partition snd gsGuesses - in go ( todo ++ cands ) giveUp ( done ++ sols ) + in go ( bis, newt + 1 ) ( todo ++ cands ) giveUp ( done ++ sols ) else -- Gauss–Seidel failed to shrink the boxes. -- Bisect along the widest dimension instead. @@ -1295,11 +1296,11 @@ intervalNewtonGS precondMethod minWidth eqs = | otherwise = [ ( t, i, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_mid ) ) , ( t, i, 𝕀 ( ℝ1 s_mid ) ( ℝ1 s_hi ) ) ] - in go ( bisGuesses ++ cands ) giveUp sols + in go ( bis + 1, newt ) ( bisGuesses ++ cands ) giveUp sols -- Box doesn't contain a solution: discard it. | otherwise - -> go cands giveUp sols + -> go ( bis, newt ) cands giveUp sols where t_mid = 0.5 * ( t_lo + t_hi ) s_mid = 0.5 * ( s_lo + s_hi ) @@ -1336,7 +1337,7 @@ precondition meth jac_mid a@( a1, a2 ) b = !a22 = 0.5 * ( a22_lo + a22_hi ) !d = a11 * a22 - a12 * a21 , not ( nearZero d ) - , let !precond = ( ℝ2 a22 -a21, ℝ2 -a12 a11 ) + , let !precond = ( ℝ2 a22 -a12, ℝ2 -a21 a11 ) !inv = recip d f x = scale inv $ matMulVec precond x -> ( ( f a1, f a2 ), f b ) @@ -1345,21 +1346,19 @@ precondition meth jac_mid a@( a1, a2 ) b = scale :: Double -> 𝕀ℝ 2 -> 𝕀ℝ 2 scale s ( 𝕀 ( ℝ2 a1_lo a2_lo ) ( ℝ2 a1_hi a2_hi ) ) - | I ( Rounded b1_lo ) ( Rounded b1_hi ) - <- I ( Rounded s ) ( Rounded s ) * I ( Rounded a1_lo ) ( Rounded a1_hi ) - , I ( Rounded b2_lo ) ( Rounded b2_hi ) - <- I ( Rounded s ) ( Rounded s ) * I ( Rounded a2_lo ) ( Rounded a2_hi ) + | 𝕀 b1_lo b1_hi <- scaleInterval s ( 𝕀 a1_lo a1_hi ) + , 𝕀 b2_lo b2_hi <- scaleInterval s ( 𝕀 a2_lo a2_hi ) = 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi ) matMulVec :: ( ℝ 2, ℝ 2 ) -> 𝕀ℝ 2 -> 𝕀ℝ 2 -matMulVec ( ℝ2 a11 a21, ℝ2 a12 a22 ) ( 𝕀 ( ℝ2 u_lo v_lo ) ( ℝ2 u_hi v_hi ) ) = - let !( I ( Rounded u'_lo ) ( Rounded u'_hi ) ) = - I ( Rounded a11 ) ( Rounded a11 ) * I ( Rounded u_lo ) ( Rounded u_hi ) - + I ( Rounded a12 ) ( Rounded a12 ) * I ( Rounded v_lo ) ( Rounded v_hi ) - !( I ( Rounded v'_lo ) ( Rounded v'_hi ) ) = - I ( Rounded a21 ) ( Rounded a21 ) * I ( Rounded u_lo ) ( Rounded u_hi ) - + I ( Rounded a22 ) ( Rounded a22 ) * I ( Rounded v_lo ) ( Rounded v_hi ) - in 𝕀 ( ℝ2 u'_lo v'_lo ) ( ℝ2 u'_hi v'_hi ) +matMulVec ( ℝ2 a11 a21, ℝ2 a12 a22 ) ( 𝕀 ( ℝ2 u_lo v_lo ) ( ℝ2 u_hi v_hi ) ) + | 𝕀 u'_lo u'_hi <- + 𝕀 a11 a11 * 𝕀 u_lo u_hi + + 𝕀 a12 a12 * 𝕀 v_lo v_hi + , 𝕀 v'_lo v'_hi <- + 𝕀 a21 a21 * 𝕀 u_lo u_hi + + 𝕀 a22 a22 * 𝕀 v_lo v_hi + = 𝕀 ( ℝ2 u'_lo v'_lo ) ( ℝ2 u'_hi v'_hi ) cmpℝ2 :: ( Double -> Double -> Bool ) -> ℝ 2 -> ℝ 2 -> Bool diff --git a/src/splines/Math/Bezier/Stroke/EnvelopeEquation.hs b/src/splines/Math/Bezier/Stroke/EnvelopeEquation.hs index 9e9c19a..41c4b3c 100644 --- a/src/splines/Math/Bezier/Stroke/EnvelopeEquation.hs +++ b/src/splines/Math/Bezier/Stroke/EnvelopeEquation.hs @@ -1,4 +1,5 @@ {-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE PolyKinds #-} {-# LANGUAGE RebindableSyntax #-} {-# LANGUAGE ScopedTypeVariables #-} @@ -11,8 +12,10 @@ module Math.Bezier.Stroke.EnvelopeEquation import Prelude hiding ( Num(..), (^) ) import Data.Coerce ( Coercible, coerce ) +import Data.Functor.Identity + ( Identity(..) ) import Data.Kind - ( Constraint ) + ( Type, Constraint ) import GHC.TypeNats ( Nat, type (-) ) @@ -25,7 +28,7 @@ import Math.Algebra.Dual import qualified Math.Bezier.Cubic as Cubic import qualified Math.Bezier.Quadratic as Quadratic import Math.Differentiable - ( type ExtentOrder ) + ( type I ) import Math.Interval import Math.Linear import Math.Module @@ -37,34 +40,35 @@ import Math.Ring -- | The value and derivative of a brush stroke at a given coordinate -- \( (t_0, s_0) \), together with the value of the envelope equation at that -- point. -data StrokeDatum i +type StrokeDatum :: Nat -> k -> Type +data StrokeDatum k i = StrokeDatum { -- | Path \( p(t_0) \). - dpath :: D ( ExtentOrder i ) ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) + dpath :: D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) -- | Brush shape \( b(t_0, s_0) \). - , dbrush :: D ( ExtentOrder i ) ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) + , dbrush :: D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) -- Everything below can be computed in terms of the first two fields. -- | Stroke \( c(t_0,s_0) = p(t_0) + b(t_0,s_0) \). - , dstroke :: D ( ExtentOrder i ) ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) + , dstroke :: D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) -- | Envelope function -- -- \[ E(t_0,s_0) = \left ( \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} \right )_{(t_0,s_0)}. \] - , ee :: D ( ExtentOrder i - 1 ) ( I i ( ℝ 2 ) ) ( I i ( ℝ 1 ) ) + , ee :: D ( k - 1 ) ( I i ( ℝ 2 ) ) ( I i ( ℝ 1 ) ) -- \[ \frac{\partial E}{\partial s} \frac{\mathrm{d} c}{\mathrm{d} t}, \] -- -- where \( \frac{\mathrm{d} c}{\mathrm{d} t} \) -- -- denotes a total derivative. - , 𝛿E𝛿sdcdt :: D ( ExtentOrder i - 2 ) ( I i ( ℝ 2 ) ) ( T ( I i ( ℝ 2 ) ) ) + , 𝛿E𝛿sdcdt :: D ( k - 2 ) ( I i ( ℝ 2 ) ) ( T ( I i ( ℝ 2 ) ) ) } -deriving stock instance Show ( StrokeDatum 'Point ) -deriving stock instance Show ( StrokeDatum 'Interval ) +deriving stock instance Show ( StrokeDatum 2 () ) +deriving stock instance Show ( StrokeDatum 3 𝕀 ) -------------------------------------------------------------------------------- @@ -123,8 +127,7 @@ class CurveOrder k where -- -- denotes a total derivative. envelopeEquation :: forall i - . ( k ~ ExtentOrder i - , D ( k - 2 ) ( I i ( ℝ 2 ) ) ~ D ( k - 2 ) ( ℝ 2 ) + . ( D ( k - 2 ) ( I i ( ℝ 2 ) ) ~ D ( k - 2 ) ( ℝ 2 ) , D ( k - 1 ) ( I i ( ℝ 2 ) ) ~ D ( k - 1 ) ( ℝ 2 ) , D k ( I i ( ℝ 2 ) ) ~ D k ( ℝ 2 ) , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) @@ -134,7 +137,6 @@ class CurveOrder k where -> ( D ( k - 1 ) ( I i ( ℝ 2 ) ) ( I i ( ℝ 1 ) ) , D ( k - 2 ) ( I i ( ℝ 2 ) ) ( T ( I i ( ℝ 2 ) ) ) ) - instance CurveOrder 2 where line ( Segment a b :: Segment b ) = D \ ( coerce -> t ) -> diff --git a/src/splines/Math/Differentiable.hs b/src/splines/Math/Differentiable.hs index 433f5b8..4f583ec 100644 --- a/src/splines/Math/Differentiable.hs +++ b/src/splines/Math/Differentiable.hs @@ -1,7 +1,8 @@ +{-# LANGUAGE PolyKinds #-} {-# LANGUAGE UndecidableInstances #-} module Math.Differentiable - ( ExtentOrder, Differentiable, DiffInterp ) + ( I, Differentiable, DiffInterp ) where -- base @@ -13,52 +14,53 @@ import GHC.TypeNats -- MetaBrush import Math.Algebra.Dual ( D, HasChainRule ) +import Math.Interval + ( 𝕀 ) import Math.Linear import Math.Module -import Math.Interval - ( Extent(..), type I ) import Math.Ring ( Transcendental ) -------------------------------------------------------------------------------- -type ExtentOrder :: Extent -> Nat -type family ExtentOrder e where - ExtentOrder 'Point = 2 - ExtentOrder 'Interval = 3 --- Currently we're doing order 2 derivatives for the brush stroke fitting, --- but order 3 derivatives for the interval Newton method to find cusps. +-- | Type family to parametrise over both pointwise and interval computations. +-- +-- Use '()' parameter for points, and '𝕀' parameter for intervals. +type I :: k -> Type -> Type +type family I i a +type instance I () a = a +type instance I 𝕀 a = 𝕀 a -type Differentiable :: Extent -> Type -> Constraint +type Differentiable :: Nat -> k -> Type -> Constraint class ( Module ( I i Double ) ( T ( I i u ) ) - , HasChainRule ( I i Double ) ( ExtentOrder i ) ( I i u ) - , Traversable ( D ( ExtentOrder i ) ( I i u ) ) - ) => Differentiable i u + , HasChainRule ( I i Double ) k ( I i u ) + , Traversable ( D k ( I i u ) ) + ) => Differentiable k i u instance ( Module ( I i Double ) ( T ( I i u ) ) - , HasChainRule ( I i Double ) ( ExtentOrder i ) ( I i u ) - , Traversable ( D ( ExtentOrder i ) ( I i u ) ) - ) => Differentiable i u + , HasChainRule ( I i Double ) k ( I i u ) + , Traversable ( D k ( I i u ) ) + ) => Differentiable k i u -type DiffInterp :: Extent -> Type -> Constraint +type DiffInterp :: Nat -> k -> Type -> Constraint class - ( Differentiable i u + ( Differentiable k i u , Interpolatable ( I i Double ) ( I i u ) , Module ( I i Double ) ( T ( I i Double ) ) - , Module ( D ( ExtentOrder i ) ( I i u ) ( I i Double ) ) - ( D ( ExtentOrder i ) ( I i u ) ( I i ( ℝ 2 ) ) ) - , Transcendental ( D ( ExtentOrder i ) ( I i u ) ( I i Double ) ) - , Applicative ( D ( ExtentOrder i ) ( I i u ) ) + , Module ( D k ( I i u ) ( I i Double ) ) + ( D k ( I i u ) ( I i ( ℝ 2 ) ) ) + , Transcendental ( D k ( I i u ) ( I i Double ) ) + , Applicative ( D k ( I i u ) ) , Representable ( I i Double ) ( I i u ) - ) => DiffInterp i u + ) => DiffInterp k i u instance - ( Differentiable i u + ( Differentiable k i u , Interpolatable ( I i Double ) ( I i u ) , Module ( I i Double ) ( T ( I i Double ) ) - , Module ( D ( ExtentOrder i ) ( I i u ) ( I i Double ) ) - ( D ( ExtentOrder i ) ( I i u ) ( I i ( ℝ 2 ) ) ) - , Transcendental ( D ( ExtentOrder i ) ( I i u ) ( I i Double ) ) - , Applicative ( D ( ExtentOrder i ) ( I i u ) ) + , Module ( D k ( I i u ) ( I i Double ) ) + ( D k ( I i u ) ( I i ( ℝ 2 ) ) ) + , Transcendental ( D k ( I i u ) ( I i Double ) ) + , Applicative ( D k ( I i u ) ) , Representable ( I i Double ) ( I i u ) - ) => DiffInterp i u + ) => DiffInterp k i u diff --git a/src/splines/Math/Interval.hs b/src/splines/Math/Interval.hs index 47f0900..7920315 100644 --- a/src/splines/Math/Interval.hs +++ b/src/splines/Math/Interval.hs @@ -4,21 +4,19 @@ {-# OPTIONS_GHC -Wno-orphans #-} +{-# OPTIONS_GHC -ddump-splices -ddump-simpl -ddump-cmm -ddump-to-file -dno-typeable-binds + -dsuppress-unfoldings -dsuppress-coercions #-} + module Math.Interval - ( 𝕀(..), 𝕀ℝ - , Extent(..), type I + ( 𝕀(𝕀), inf, sup + , scaleInterval + , 𝕀ℝ , singleton, nonDecreasing ) where -- base import Prelude hiding ( Num(..) ) -import Data.Coerce - ( coerce ) -import Data.Kind - ( Type ) -import Data.Monoid - ( Sum(..) ) -- acts import Data.Act @@ -37,7 +35,7 @@ import Math.Algebra.Dual import Math.Algebra.Dual.Internal ( chainRule1NQ ) import Math.Interval.Internal - ( 𝕀(..) ) + ( 𝕀(𝕀), inf, sup, scaleInterval ) import Math.Linear ( ℝ(..), T(..) , RepresentableQ(..) @@ -52,14 +50,6 @@ import Math.Ring type 𝕀ℝ n = 𝕀 ( ℝ n ) type instance D k ( 𝕀 v ) = D k v --- Handling points and intervals uniformly. -data Extent = Point | Interval - -type I :: Extent -> Type -> Type -type family I i a where - I 'Point a = a - I 'Interval a = 𝕀 a - singleton :: a -> 𝕀 a singleton a = 𝕀 a a @@ -76,9 +66,9 @@ deriving via ViaAbelianGroup ( T ( 𝕀 Double ) ) instance Group ( T ( 𝕀 Double ) ) instance Act ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where - T g • a = coerce ( Sum g • a ) + T g • a = g + a instance Torsor ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where - a --> b = T $ getSum ( a --> b ) + a --> b = T $ b - a ------------------------------------------------------------------------------- @@ -171,9 +161,9 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 1 ) where | Just i <- isLinear mon -> [|| f $$( tabulateQ \ j -> if | j == i - -> [|| 1 ||] + -> [|| 𝕀 1 1 ||] | otherwise - -> [|| 0 ||] + -> [|| 𝕀 0 0 ||] ) ||] | otherwise -> [|| unT o ||] @@ -206,9 +196,9 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 1 ) where | Just i <- isLinear mon -> [|| f $$( tabulateQ \ j -> if | j == i - -> [|| 1 ||] + -> [|| 𝕀 1 1 ||] | otherwise - -> [|| 0 ||] + -> [|| 𝕀 0 0 ||] ) ||] | otherwise -> [|| unT o ||] @@ -241,9 +231,9 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 2 ) where | Just i <- isLinear mon -> [|| f $$( tabulateQ \ j -> if | j == i - -> [|| 1 ||] + -> [|| 𝕀 1 1 ||] | otherwise - -> [|| 0 ||] + -> [|| 𝕀 0 0 ||] ) ||] | otherwise -> [|| unT o ||] @@ -276,9 +266,9 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 2 ) where | Just i <- isLinear mon -> [|| f $$( tabulateQ \ j -> if | j == i - -> [|| 1 ||] + -> [|| 𝕀 1 1 ||] | otherwise - -> [|| 0 ||] + -> [|| 𝕀 0 0 ||] ) ||] | otherwise -> [|| unT o ||] @@ -311,9 +301,9 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 3 ) where | Just i <- isLinear mon -> [|| f $$( tabulateQ \ j -> if | j == i - -> [|| 1 ||] + -> [|| 𝕀 1 1 ||] | otherwise - -> [|| 0 ||] + -> [|| 𝕀 0 0 ||] ) ||] | otherwise -> [|| unT o ||] @@ -346,9 +336,9 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 3 ) where | Just i <- isLinear mon -> [|| f $$( tabulateQ \ j -> if | j == i - -> [|| 1 ||] + -> [|| 𝕀 1 1 ||] | otherwise - -> [|| 0 ||] + -> [|| 𝕀 0 0 ||] ) ||] | otherwise -> [|| unT o ||] @@ -381,9 +371,9 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 4 ) where | Just i <- isLinear mon -> [|| f $$( tabulateQ \ j -> if | j == i - -> [|| 1 ||] + -> [|| 𝕀 1 1 ||] | otherwise - -> [|| 0 ||] + -> [|| 𝕀 0 0 ||] ) ||] | otherwise -> [|| unT o ||] @@ -416,9 +406,9 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 4 ) where | Just i <- isLinear mon -> [|| f $$( tabulateQ \ j -> if | j == i - -> [|| 1 ||] + -> [|| 𝕀 1 1 ||] | otherwise - -> [|| 0 ||] + -> [|| 𝕀 0 0 ||] ) ||] | otherwise -> [|| unT o ||] diff --git a/src/splines/Math/Interval/FMA.hs b/src/splines/Math/Interval/FMA.hs new file mode 100644 index 0000000..eea05f6 --- /dev/null +++ b/src/splines/Math/Interval/FMA.hs @@ -0,0 +1,152 @@ +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} + +module Math.Interval.FMA ( addI, subI, prodI, divI, posPowI ) where + +-- base +import Data.Bits + ( Bits((.&.), shiftL) ) +import Data.Word + ( Word32, Word64 ) +import GHC.Exts + ( Double(D#), fmsubDouble#, fnmaddDouble# ) +import GHC.Float + ( castFloatToWord32 , castWord32ToFloat + , castDoubleToWord64, castWord64ToDouble + ) + +-------------------------------------------------------------------------------- + +class ( RealFloat f, Num b, Bits b ) => FPBits f b | f -> b, b -> f where + toBits :: f -> b + fromBits :: b -> f + -- | Size in bytes. + sizeOf :: Int + +instance FPBits Float Word32 where + toBits = castFloatToWord32 + fromBits = castWord32ToFloat + sizeOf = 4 + +instance FPBits Double Word64 where + toBits = castDoubleToWord64 + fromBits = castWord64ToDouble + sizeOf = 8 + + +{-# SPECIALISE nextAfter :: Float -> Float -> Float #-} +{-# SPECIALISE nextAfter :: Double -> Double -> Double #-} +{-# INLINEABLE nextAfter #-} +-- | @nextAfter a b@ computes the next floating-point value after @a@ +-- in the direction of @b@. +nextAfter :: forall f b. FPBits f b => f -> f -> f +nextAfter a b + | isNaN a + = a + | isNaN b + = b + | a == b + = b + | otherwise + = let !res_bits + | a == 0 + , let !sgn_mask = 1 `shiftL` ( sizeOf @f * 8 - 1 ) + = ( toBits b .&. sgn_mask ) + 1 + | ( a < b ) == ( a > 0 ) + = toBits a + 1 + | otherwise + = toBits a - 1 + in fromBits res_bits + +{-# SPECIALISE succFP :: Float -> Float #-} +{-# SPECIALISE succFP :: Double -> Double #-} +{-# INLINEABLE succFP #-} +-- | The next floating-point number. +succFP :: forall f b. FPBits f b => f -> f +succFP x = nextAfter x (1/0) + + +{-# SPECIALISE prevFP :: Float -> Float #-} +{-# SPECIALISE prevFP :: Double -> Double #-} +{-# INLINEABLE prevFP #-} +-- | The previous floating-point number. +prevFP :: forall f b. FPBits f b => f -> f +prevFP x = nextAfter x (-1/0) + +-------------------------------------------------------------------------------- + +{-# INLINE withError #-} +withError :: Double -> Double -> ( Double, Double ) +withError f e = + case compare e 0 of + LT -> ( prevFP f, f ) + EQ -> ( f, f ) + GT -> ( f, succFP f ) + +addI :: Double -> Double -> ( Double, Double ) +addI a b = + let !s = a + b + !a' = s - b + !b' = s - a' + !da = a - a' + !db = b - b' + !e = da + db + in s `withError` e + +subI :: Double -> Double -> ( Double, Double ) +subI a b = + let !s = a - b + !a' = s - b + !b' = s - a' + !da = a - a' + !db = b - b' + !e = da + db + in s `withError` e + +{-# INLINE fmsubDouble #-} +fmsubDouble :: Double -> Double -> Double -> Double +fmsubDouble ( D# x ) ( D# y ) ( D# z ) = D# ( fmsubDouble# x y z ) +{-# INLINE fnmaddDouble #-} +fnmaddDouble :: Double -> Double -> Double -> Double +fnmaddDouble ( D# x ) ( D# y ) ( D# z ) = D# ( fnmaddDouble# x y z ) + +prodI :: Double -> Double -> ( Double, Double ) +prodI a b = + let !p = a * b + !e = fmsubDouble a b p + in p `withError` e + +divI :: Double -> Double -> ( Double, Double ) +divI a b = + let !r = a / b + !e = fnmaddDouble r b a + in r `withError` e + +-- | Power of a **non-negative** number to a natural power. +posPowI :: Double -- ^ Assumed to be non-negative! + -> Word -> ( Double, Double ) +posPowI _ 0 = ( 1, 1 ) +posPowI f 1 = ( f, f ) +posPowI f 2 = prodI f f +posPowI f n + | even n + , let m = n `quot` 2 + ( f²_lo, f²_hi ) = prodI f f + = ( fst $ posPowI f²_lo m, snd $ posPowI f²_hi m ) + | otherwise + , let m = n `quot` 2 + ( f²_lo, f²_hi ) = prodI f f + = ( fst $ posPowAcc f²_lo m f, snd $ posPowAcc f²_hi m f ) + +posPowAcc :: Double -> Word -> Double -> ( Double, Double ) +posPowAcc f 1 x = prodI f x +posPowAcc f n x + | even n + , let m = n `quot` 2 + ( f²_lo, f²_hi ) = prodI f f + = ( fst $ posPowAcc f²_lo m x, snd $ posPowAcc f²_hi m x ) + | otherwise + , let m = n `quot` 2 + ( f²_lo, f²_hi ) = prodI f f + ( y_lo, y_hi ) = prodI f x + = ( fst $ posPowAcc f²_lo m y_lo, snd $ posPowAcc f²_hi m y_hi ) diff --git a/src/splines/Math/Interval/Internal.hs b/src/splines/Math/Interval/Internal.hs index c70ebcb..a8d692a 100644 --- a/src/splines/Math/Interval/Internal.hs +++ b/src/splines/Math/Interval/Internal.hs @@ -1,11 +1,12 @@ -{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE UndecidableInstances #-} {-# OPTIONS_GHC -Wno-orphans #-} module Math.Interval.Internal - ( 𝕀(.., 𝕀) ) + ( 𝕀(𝕀), inf, sup, scaleInterval ) where -- base @@ -18,17 +19,24 @@ import GHC.Show -- deepseq import Control.DeepSeq - ( NFData ) + ( NFData(..) ) -- rounded-hw import Numeric.Rounded.Hardware ( Rounded(..) ) -import Numeric.Rounded.Hardware.Interval.NonEmpty - ( Interval ) import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval ( Interval(..) ) + +#ifdef USE_FMA +-- MetaBrush +import Math.Interval.FMA + ( addI, subI, prodI, divI, posPowI ) + +#else +-- rounded-hw import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval - ( sup, inf, powInt ) + ( powInt ) +#endif -- MetaBrush import Math.Linear @@ -42,29 +50,104 @@ import Math.Ring -------------------------------------------------------------------------------- -- Intervals. -newtype 𝕀 a = MkI { ival :: Interval a } +#ifdef USE_FMA +data 𝕀 a = 𝕀 { inf, sup :: !a } + +instance NFData a => NFData ( 𝕀 a ) where + rnf ( 𝕀 lo hi ) = rnf lo `seq` rnf hi + +instance Prelude.Num ( 𝕀 Double ) where + 𝕀 x_lo x_hi + 𝕀 y_lo y_hi + | let !z_lo = fst $ addI x_lo y_lo + !z_hi = snd $ addI x_hi y_hi + = 𝕀 z_lo z_hi + 𝕀 x_lo x_hi - 𝕀 y_lo y_hi + | let !z_lo = fst $ subI x_lo y_hi + !z_hi = snd $ subI x_hi y_lo + = 𝕀 z_lo z_hi + negate (𝕀 lo hi) = 𝕀 (Prelude.negate hi) (Prelude.negate lo) + (*) = (*) + fromInteger i = + let !j = Prelude.fromInteger i + in 𝕀 j j + abs (𝕀 lo hi) + | 0 <= lo + = 𝕀 lo hi + | hi <= 0 + = 𝕀 (Prelude.negate hi) (Prelude.negate lo) + | otherwise + = 𝕀 0 (max (Prelude.negate lo) hi) + signum _ = error "No implementation of signum for intervals" + +instance Ring ( 𝕀 Double ) where + 𝕀 lo1 hi1 * 𝕀 lo2 hi2 + | let !( x_min, x_max ) = prodI lo1 lo2 + !( y_min, y_max ) = prodI lo1 hi2 + !( z_min, z_max ) = prodI hi1 lo2 + !( w_min, w_max ) = prodI hi1 hi2 + = 𝕀 ( min ( min x_min y_min ) ( min z_min w_min ) ) + ( max ( max x_max y_max ) ( max z_max w_max ) ) + _ ^ 0 = 𝕀 1 1 + iv ^ 1 = iv + 𝕀 lo hi ^ n + | odd n || 0 <= lo + , let !lo' = fst $ posPowI lo n + !hi' = snd $ posPowI hi n + = 𝕀 lo' hi' + | hi <= 0 + , let !lo' = fst $ posPowI (negate hi) n + !hi' = snd $ posPowI (negate lo) n + = 𝕀 lo' hi' + | otherwise + , let !hi1 = snd $ posPowI (negate lo) n + , let !hi2 = snd $ posPowI hi n + = 𝕀 0 ( max hi1 hi2 ) + +instance Prelude.Fractional ( 𝕀 Double ) where + fromRational r = + let q = Prelude.fromRational r + in 𝕀 q q + recip (𝕀 lo hi) +-- #ifdef ASSERTS + | lo >= 0 || hi <= 0 +-- #endif + = 𝕀 ( fst $ divI 1 hi ) ( snd $ divI 1 lo ) +-- #ifdef ASSERTS + | otherwise + = error "BAD interval recip; should use extendedRecip instead" +-- #endif + _ / _ = error "TODO: interval division is not implemented" + +instance Transcendental ( 𝕀 Double ) where + pi = 𝕀 3.141592653589793 3.1415926535897936 + cos = withHW Prelude.cos + sin = withHW Prelude.sin + +{-# INLINE withHW #-} +withHW :: (Interval.Interval a -> Interval.Interval b) -> 𝕀 a -> 𝕀 b +withHW f = \ ( 𝕀 lo hi ) -> + case f ( Interval.I ( Rounded lo ) ( Rounded hi ) ) of + Interval.I ( Rounded x ) ( Rounded y ) -> 𝕀 x y + +scaleInterval :: Double -> 𝕀 Double -> 𝕀 Double +scaleInterval s ( 𝕀 lo hi ) = + case compare s 0 of + LT -> 𝕀 ( fst $ prodI s hi ) ( snd $ prodI s lo ) + EQ -> 𝕀 0 0 + GT -> 𝕀 ( fst $ prodI s lo ) ( snd $ prodI s hi ) +#else + +newtype 𝕀 a = MkI { ival :: Interval.Interval a } deriving newtype ( Prelude.Num, Prelude.Fractional, Prelude.Floating ) deriving newtype NFData -instance Eq a => Eq ( 𝕀 a ) where - 𝕀 a b == 𝕀 c d = - a == c && b == d - {-# COMPLETE 𝕀 #-} pattern 𝕀 :: a -> a -> 𝕀 a pattern 𝕀 x y = MkI ( Interval.I ( Rounded x ) ( Rounded y ) ) -instance Show a => Show ( 𝕀 a ) where - showsPrec _ ( 𝕀 x y ) - = showString "[" - . showsPrec 0 x - . showCommaSpace - . showsPrec 0 y - . showString "]" -deriving via ViaPrelude ( 𝕀 Double ) - instance AbelianGroup ( 𝕀 Double ) -deriving via ViaPrelude ( 𝕀 Double ) - instance AbelianGroup ( T ( 𝕀 Double ) ) +inf, sup :: 𝕀 a -> a +inf (𝕀 a _) = a +sup (𝕀 _ b) = b instance Ring ( 𝕀 Double ) where MkI i1 * MkI i2 = MkI $ i1 Prelude.* i2 @@ -74,16 +157,40 @@ instance Ring ( 𝕀 Double ) where -- accidentally use (^) from Prelude. deriving via ViaPrelude ( 𝕀 Double ) - instance Field ( 𝕀 Double ) + instance Transcendental ( 𝕀 Double ) + +scaleInterval :: Double -> 𝕀 Double -> 𝕀 Double +scaleInterval s iv = 𝕀 s s * iv -- TODO: could be better +#endif deriving via ViaPrelude ( 𝕀 Double ) - instance Transcendental ( 𝕀 Double ) + instance AbelianGroup ( 𝕀 Double ) +deriving via ViaPrelude ( 𝕀 Double ) + instance AbelianGroup ( T ( 𝕀 Double ) ) +deriving via ViaPrelude ( 𝕀 Double ) + instance Field ( 𝕀 Double ) + +-------------------------------------------------------------------------------- + +instance Eq a => Eq ( 𝕀 a ) where + 𝕀 a b == 𝕀 c d = + a == c && b == d + +instance Show a => Show ( 𝕀 a ) where + showsPrec _ ( 𝕀 x y ) + = showString "[" + . showsPrec 0 x + . showCommaSpace + . showsPrec 0 y + . showString "]" + +-------------------------------------------------------------------------------- type instance RepDim ( 𝕀 u ) = RepDim u instance RepresentableQ r u => RepresentableQ ( 𝕀 r ) ( 𝕀 u ) where tabulateQ f = - let !lo = tabulateQ @r @u ( \ i -> [|| getRounded $ Interval.inf $ ival $$( f i ) ||] ) - !hi = tabulateQ @r @u ( \ i -> [|| getRounded $ Interval.sup $ ival $$( f i ) ||] ) + let !lo = tabulateQ @r @u ( \ i -> [|| inf $ $$( f i ) ||] ) + !hi = tabulateQ @r @u ( \ i -> [|| sup $ $$( f i ) ||] ) in [|| 𝕀 $$lo $$hi ||] indexQ d i = @@ -95,8 +202,8 @@ instance RepresentableQ r u => RepresentableQ ( 𝕀 r ) ( 𝕀 u ) where ||] instance Representable r u => Representable ( 𝕀 r ) ( 𝕀 u ) where tabulate f = - let !lo = tabulate @r @u ( \ i -> getRounded $ Interval.inf $ ival ( f i ) ) - !hi = tabulate @r @u ( \ i -> getRounded $ Interval.sup $ ival ( f i ) ) + let !lo = tabulate @r @u ( \ i -> inf $ f i ) + !hi = tabulate @r @u ( \ i -> sup $ f i ) in 𝕀 lo hi index d i =