improve rotations in interval arithmetic

Computing a rotation in interval arithmetic can lose tightness.
Instead of computing

   a cos phi + b sin phi

which doesn't account for the difference in phase in the two sinusoids,
it is better to use

  sqrt (a² + b²) * cos ( phi - atan2(a,b) )

which correctly estimates the maximum amplitude of the sum to be
sqrt(a²+b²) instead of abs(a) + abs(b).

This seems to worsen the performance of the interval Newton method
at the moment, possibly due to the complexity of the new formula,
which involves computing atan(b/a). Further investigation is needed.
This commit is contained in:
sheaf 2023-05-14 21:38:25 +02:00
parent 50d99e1e4b
commit 96aa38b3c3
12 changed files with 796 additions and 212 deletions

View file

@ -273,6 +273,9 @@ executable cusps
hs-source-dirs:
src/cusps
default-language:
Haskell2010
main-is:
Main.hs

View file

@ -167,9 +167,9 @@ runApplication application = do
-- ]
-- }
Spline
{ splineStart = mkPoint ( 2 10 -20 ) 2 1 0
{ splineStart = mkPoint ( 2 0 0 ) 1 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 100 0 ) 10 10 (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 }
]

View file

@ -1,6 +1,7 @@
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE PartialTypeSignatures #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE RebindableSyntax #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE UndecidableInstances #-}
@ -10,7 +11,11 @@
module Main where
-- base
import Prelude hiding ( Num(..), Fractional(..), Floating(..), (^) )
import Prelude
hiding
( Num(..), Fractional(..), Floating(..), (^), (/)
, fromInteger, fromRational
)
import qualified Prelude
import GHC.Generics
( Generic )
@ -29,21 +34,29 @@ import Data.Generics.Product.Typed
-- MetaBrush
import Math.Algebra.Dual
import qualified Math.Bezier.Quadratic as Quadratic
import qualified Math.Bezier.Cubic as Cubic
import Math.Bezier.Spline
import Math.Bezier.Stroke
( brushStrokeData, pathAndUsedParams )
import Math.Bezier.Stroke.EnvelopeEquation
( StrokeDatum(..) )
( HasBézier(..), StrokeDatum(..) )
import Math.Differentiable
( I )
import Math.Linear
( (..), T(..)
( (..), T(..), Segment(..)
, Fin(..)
, Representable(..) )
, Representable(..)
)
import Math.Module
( Module(..), Cross(..) )
( Module(..), Cross(..)
, lerp
)
import Math.Ring
( AbelianGroup(..), Ring(..), Transcendental(..) )
( AbelianGroup(..), Ring(..)
, Field(..), Floating(..), Transcendental(..)
, ifThenElse
)
import Math.Interval.Abstract
@ -66,16 +79,19 @@ outlineFunction
, 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 ) )
, Module Double (T brushParams)
, Torsor (T brushParams) brushParams
, Representable Double 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 ) )
, HasBézier 3 i
)
=> ( I i Double -> I i ( 1 ) )
-> ( I i ( 1 ) -> I i Double )
@ -97,11 +113,15 @@ outlineFunction co1 co2 single brush brush_index sp0 crv = strokeData
strokeData :: I i ( 1 ) -> I i ( 1 ) -> StrokeDatum 3 i
strokeData =
fmap ( `Seq.index` brush_index ) $
brushStrokeData @3 @brushParams @i co1 co2 path params brush
brushStrokeData @3 @brushParams co1 co2 path params brush
boo :: ( HasType ( 2 ) ( PointData brushParams )
, Show brushParams
, Representable Double brushParams
, Module Double ( T brushParams )
, Module ( AI Double ) ( T ( AI brushParams ) )
, Torsor ( T brushParams ) brushParams
, Torsor ( T ( AI brushParams ) ) ( AI brushParams )
, HasChainRule ( AI Double ) 3 ( AI brushParams )
, Traversable ( D 3 brushParams )
@ -118,9 +138,9 @@ boo x i y z
ellipseTest :: StrokeDatum 3 AI
ellipseTest =
boo ( ellipseBrush @AI ( Pt . Val ) scaleI ) 1
( mkPoint ( 2 10 -20 ) 2 1 0 )
( LineTo { curveEnd = NextPoint $ mkPoint ( 2 10 10 ) 10 5 ( pi Prelude./ 4 )
boo ( ellipseBrush ( Pt . Val ) scaleI ) 1
( mkPoint ( 2 0 0 ) 1.0111 1.0222 0 )
( LineTo { curveEnd = NextPoint $ mkPoint ( 2 100 0 ) 10.0333 10.0444 ( pi Prelude./ 2 )
, curveData = () } )
where
mkPoint :: 2 -> Double -> Double -> Double -> PointData ( 3 )
@ -214,6 +234,28 @@ instance HasChainRule ( AI Double ) 3 ( AI ( 4 ) ) where
[|| o ||] [|| p ||] [|| s ||]
[|| df ||] [|| dg ||] )
instance HasBézier 3 AI where
line co ( Segment a b :: Segment b ) =
D \ ( co -> t ) ->
D31 ( lerp @( T b ) t a b )
( a --> b )
origin
origin
bezier2 co ( bez :: Quadratic.Bezier b ) =
D \ ( co -> t ) ->
D31 ( EvalBez2 bez t )
( Quadratic.bezier' bez t )
( Quadratic.bezier'' bez )
origin
bezier3 co ( bez :: Cubic.Bezier b ) =
D \ ( co -> t ) ->
D31 ( EvalBez3 bez t )
( T $ EvalBez2 ( fmap unT $ Cubic.derivative ( fmap T bez ) ) t )
( Cubic.bezier'' bez t )
( Cubic.bezier''' bez )
--------------------------------------------------------------------------------
-- Brushes (TODO copied from MetaBrush.Asset.Brushes)
@ -237,40 +279,60 @@ circleSpline p = sequenceA $
Bezier3To ( p κ -1 ) ( p 1 -κ ) BackToStart ()
{-# INLINE circleSpline #-}
ellipseBrush :: forall {t} (i :: t) k irec
. ( irec ~ I i ( 3 )
ellipseBrush :: forall k irec
. ( irec ~ AI ( 3 )
, Module
( D k irec ( I i Double ) )
( D k irec ( I i ( 2 ) ) )
, Module ( I i Double ) ( T ( I i Double ) )
, HasChainRule ( I i Double ) k irec
, Representable ( I i Double ) irec
( D k irec ( AI Double ) )
( D k irec ( AI ( 2 ) ) )
, Module ( AI Double ) ( T ( AI Double ) )
, HasChainRule ( AI Double ) k irec
, Representable ( AI Double ) irec
, Applicative ( D k irec )
, Transcendental ( D k irec ( I i Double ) )
, Transcendental ( D k irec ( AI Double ) )
)
=> ( forall a. a -> I i a )
-> ( Double -> I i Double -> I i Double )
-> C k irec ( Spline 'Closed () ( I i ( 2 ) ) )
=> ( forall a. a -> AI a )
-> ( Double -> AI Double -> AI Double )
-> C k irec ( Spline 'Closed () ( AI ( 2 ) ) )
ellipseBrush mkI scaleI =
D \ params ->
let a, b, phi :: D k irec ( I i Double )
let a, b, phi :: D k irec ( AI Double )
a = runD ( var @_ @k ( Fin 1 ) ) params
b = runD ( var @_ @k ( Fin 2 ) ) params
phi = runD ( var @_ @k ( Fin 3 ) ) params
mkPt :: Double -> Double -> D k irec ( I i ( 2 ) )
mkPt :: Double -> Double -> D k irec ( AI ( 2 ) )
mkPt x y
= let !a' = mul x a
!b' = mul y b
= let !x' = a `scaledBy` x
!y' = b `scaledBy` y
{-
!c = cos phi
!s = sin phi
in ( a' * c - b' * s ) *^ e_x
^+^ ( a' * s + b' * c ) *^ e_y
in circleSpline @i @k @( 3 ) @( 2 ) mkPt
in ( x' * c - y' * s ) *^ e_x
^+^ ( x' * s + y' * c ) *^ e_y
-}
-- {-
!r = sqrt $ x' ^ 2 + y' ^ 2
!arctgt = atan ( y' / x' )
-- a and b are always strictly positive, so we can compute
-- the quadrant using only x and y, which are constants.
!theta
| x > 0
= arctgt
| x < 0
= if y >= 0 then arctgt + pi else arctgt - pi
| otherwise
= if y >= 0 then 0.5 * pi else -0.5 * pi
!phi' = phi + theta
in
( r * cos phi' ) *^ e_x
^+^ ( r * sin phi' ) *^ e_y
-- -}
in circleSpline @AI @k @( 3 ) @( 2 ) mkPt
where
e_x, e_y :: D k irec ( I i ( 2 ) )
e_x, e_y :: D k irec ( AI ( 2 ) )
e_x = pure $ mkI $ 2 1 0
e_y = pure $ mkI $ 2 0 1
mul :: Double -> D k irec ( I i Double ) -> D k irec ( I i Double )
mul k = fmap ( scaleI k )
scaledBy :: D k irec ( AI Double ) -> Double -> D k irec ( AI Double )
scaledBy d k = fmap ( scaleI k ) d
{-# INLINEABLE ellipseBrush #-}

View file

@ -13,13 +13,21 @@ import Prelude hiding
import qualified Prelude
( Num(..), Fractional(..) )
import Data.Coerce
( Coercible, coerce )
( coerce )
import Data.Dynamic
( Dynamic, toDyn, fromDynamic, dynTypeRep )
import Data.Foldable
( toList )
import Data.List
( intersperse )
import Data.Proxy
( Proxy(..) )
import Data.Typeable
( Typeable, typeRep )
import GHC.Show
( showCommaSpace )
( showSpace, showCommaSpace )
import GHC.TypeNats
( KnownNat )
-- acts
import Data.Act
@ -28,8 +36,8 @@ import Data.Act
-- containers
import Data.Map
( Map )
import qualified Data.Map as Map
( (!) )
import qualified Data.Map.Strict as Map
( empty, insertLookupWithKey, lookup )
-- groups
import Data.Group
@ -37,6 +45,8 @@ import Data.Group
-- MetaBrush
import Math.Algebra.Dual
import qualified Math.Bezier.Quadratic as Quadratic
import qualified Math.Bezier.Cubic as Cubic
import Math.Bezier.Stroke.EnvelopeEquation
( StrokeDatum(..) )
import Math.Differentiable
@ -50,21 +60,24 @@ import Math.Linear
import Math.Module
( Module(..), Cross(..) )
import Math.Ring
( AbelianGroup(..), Ring(..), Field(..), Transcendental(..)
( AbelianGroup(..), Ring(..), Field(..)
, Floating(..), Transcendental(..)
, ViaPrelude(..)
)
--------------------------------------------------------------------------------
-- EDSL for inspection
type instance D k (AI a) = D k (𝕀 a)
type instance D k ( AI a ) = D k ( 𝕀 a )
-- TODO: I think we want a custom representation for D k (AI a),
-- so that we can recover sharing.
type instance I AI a = AI a
-- | A value: a constant, or a variable name.
data Value a where
Val :: a -> Value a
Var :: String -> Value Double
Var :: ( Typeable a, Eq a, Show a ) => String -> Value a
-- | A scalar or a vector.
data SV a where
@ -87,16 +100,20 @@ data AI a where
Recip :: AI Double -> AI Double
Scale :: IsVec v => AI Double -> AI v -> AI v
Cross :: AI ( 2 ) -> AI ( 2 ) -> AI Double
Sqrt :: AI Double -> AI Double
Cos :: AI Double -> AI Double
Sin :: AI Double -> AI Double
Atan :: AI Double -> AI Double
Pi :: AI Double
Coerce1 :: AI Double -> AI ( 1 )
Coerce2 :: AI ( 1 ) -> AI Double
Tabulate :: ( Eq ( n ), Show ( n ), Representable Double ( n ) ) => Vec n ( AI Double ) -> AI ( n )
Index :: ( Eq ( n ), Show ( n ), Representable Double ( n ) ) => AI ( n ) -> Fin n -> AI Double
EvalBez2 :: Quadratic.Bezier (AI a) -> AI Double -> AI a
EvalBez3 :: Cubic.Bezier (AI a) -> AI Double -> AI a
Tabulate :: ( KnownNat n, Eq ( n ), Show ( n ), Representable Double ( n ) ) => Vec n ( AI Double ) -> AI ( n )
Index :: ( KnownNat n, Eq ( n ), Show ( n ), Representable Double ( n ) ) => AI ( n ) -> Fin n -> AI Double
instance ( Show a, Eq a ) => Show ( AI a ) where
showsPrec prec ai = case fst $ normalise ai of
instance ( Show a, Eq a, Typeable a ) => Show ( AI a ) where
showsPrec prec ai = case fst $ normalise Map.empty ai of
Pt v -> showsPrec prec v
Ival lo hi -> showString "[" . showsPrec 0 lo . showCommaSpace . showsPrec 0 hi . showString "]"
(:+:) Scalar iv1 iv2 -> showParen (prec > 6) $ showsPrec 6 iv1 . showString " + " . showsPrec 7 iv2
@ -113,8 +130,12 @@ instance ( Show a, Eq a ) => Show ( AI a ) where
Recip iv -> showParen (prec > 10) $ showString " 1 / " . showsPrec 11 iv
Coerce1 iv -> showsPrec prec iv
Coerce2 iv -> showsPrec prec iv
EvalBez2 bez t -> showParen (prec > 10) $ showString "bez2" . showSpace . showsPrec 11 t . showSpace . showString "[" . showBez2 bez . showString "]"
EvalBez3 bez t -> showParen (prec > 10) $ showString "bez3" . showSpace . showsPrec 11 t . showSpace . showString "[" . showBez3 bez . showString "]"
Sqrt iv -> showParen (prec > 10) $ showString "sqrt " . showsPrec 11 iv
Cos iv -> showParen (prec > 10) $ showString "cos " . showsPrec 11 iv
Sin iv -> showParen (prec > 10) $ showString "sin " . showsPrec 11 iv
Atan iv -> showParen (prec > 10) $ showString "atan " . showsPrec 11 iv
Pi -> showString "pi"
Tabulate vs -> showString "(" . foldr (.) id ( intersperse showCommaSpace fields ) . showString ")"
where
@ -123,6 +144,9 @@ instance ( Show a, Eq a ) => Show ( AI a ) where
| v <- toList vs
]
showBez2 (Quadratic.Bezier p0 p1 p2 ) = foldr (.) id ( intersperse showCommaSpace $ map (showsPrec 0) [p0, p1, p2] )
showBez3 ( Cubic.Bezier p0 p1 p2 p3) = foldr (.) id ( intersperse showCommaSpace $ map (showsPrec 0) [p0, p1, p2, p3] )
infixl 6 :+:
infixl 6 :-:
infixl 7 :*:
@ -153,15 +177,29 @@ index_maybe expr i = case expr of
-> Just $ v ! i
Coerce1 v
-> Just $ v
EvalBez2 bez t
| Just bez_i <- traverse ( `index_maybe` i ) bez
-> Just $ EvalBez2 bez_i t
EvalBez3 bez t
| Just bez_i <- traverse ( `index_maybe` i ) bez
-> Just $ EvalBez3 bez_i t
_ -> Nothing
ival :: 𝕀 a -> AI a
ival ( 𝕀 lo hi ) = Ival ( Val lo ) ( Val hi )
normalise :: Eq a => AI a -> ( AI a, Maybe ( 𝕀 a ) )
normalise expr = case expr of
Pt ( Val v ) -> ( expr, Just $ 𝕀 v v )
Pt {} -> ( expr, Nothing )
normalise :: forall a. ( Eq a, Typeable a )
=> Map String Dynamic -> AI a -> ( AI a, Maybe ( 𝕀 a ) )
normalise vars expr = case expr of
Pt ( Val v )
-> ( expr, Just $ 𝕀 v v )
Pt ( Var v )
| Just d <- Map.lookup v vars
, Just a <- fromDynamic @( AI a ) d
, Just r <- snd $ normalise vars a
-> ( ival r, Just r )
| otherwise
-> ( expr, Nothing )
Ival ( Val v1 ) ( Val v2 )
| v1 == v2
-> ( Pt ( Val v1 ), Just $ 𝕀 v1 v2 )
@ -191,16 +229,16 @@ normalise expr = case expr of
-> ( expr', Nothing )
where
expr' = (:+:) sv iv1' iv2'
(iv1', mv1) = normalise iv1
(iv2', mv2) = normalise iv2
(iv1', mv1) = normalise vars iv1
(iv2', mv2) = normalise vars iv2
(:-:) sv iv1 iv2
| Scalar <- sv
, Just ( 𝕀 0 0 ) <- mv1
-> normalise $ Negate sv iv2
-> normalise vars $ Negate sv iv2
| Vector <- sv
, Just o <- mv1
, T o == origin
-> normalise $ Negate sv iv2
-> normalise vars $ Negate sv iv2
| Scalar <- sv
, Just ( 𝕀 0 0 ) <- mv2
-> ( iv1', mv1 )
@ -216,8 +254,8 @@ normalise expr = case expr of
-> ( expr', Nothing )
where
expr' = (:-:) sv iv1' iv2'
(iv1', mv1) = normalise iv1
(iv2', mv2) = normalise iv2
(iv1', mv1) = normalise vars iv1
(iv2', mv2) = normalise vars iv2
Negate sv iv
| Just v <- mv
-> let r = case sv of { Scalar -> unT $ negate $ T v; Vector -> unT $ origin ^-^ T v }
@ -226,10 +264,14 @@ normalise expr = case expr of
-> ( expr', Nothing )
where
expr' = Negate sv iv'
(iv', mv) = normalise iv
(iv', mv) = normalise vars iv
iv1 :*: iv2
| Just ( 𝕀 0 0 ) <- mv1
-> zero
| Just ( 𝕀 1 1 ) <- mv1
-> (iv2', mv2)
| Just ( 𝕀 1 1 ) <- mv2
-> (iv1', mv1)
| Just ( 𝕀 0 0 ) <- mv2
-> zero
| Just v1 <- mv1
@ -241,8 +283,8 @@ normalise expr = case expr of
where
zero = ( ival $ 𝕀 0 0, Just $ 𝕀 0 0 )
expr' = (:*:) iv1' iv2'
(iv1', mv1) = normalise iv1
(iv2', mv2) = normalise iv2
(iv1', mv1) = normalise vars iv1
(iv2', mv2) = normalise vars iv2
iv :^: n
| Just v <- mv
-> let r = v ^ n
@ -251,7 +293,7 @@ normalise expr = case expr of
-> ( expr', Nothing )
where
expr' = iv' :^: n
(iv', mv) = normalise iv
(iv', mv) = normalise vars iv
Recip iv
| Just v <- mv
-> let r = recip v
@ -260,7 +302,7 @@ normalise expr = case expr of
-> ( expr', Nothing )
where
expr' = Recip iv'
(iv', mv) = normalise iv
(iv', mv) = normalise vars iv
Scale k iv
| Just ( 𝕀 0 0 ) <- mk
-> ( ival $ unT $ origin, Just $ unT origin )
@ -277,8 +319,8 @@ normalise expr = case expr of
-> ( expr', Nothing )
where
expr' = Scale k' iv'
(k' , mk) = normalise k
(iv', mv) = normalise iv
(k' , mk) = normalise vars k
(iv', mv) = normalise vars iv
Cross iv1 iv2
| Just v1 <- mv1
, Just v2 <- mv2
@ -288,8 +330,17 @@ normalise expr = case expr of
-> ( expr', Nothing )
where
expr' = Cross iv1' iv2'
(iv1', mv1) = normalise iv1
(iv2', mv2) = normalise iv2
(iv1', mv1) = normalise vars iv1
(iv2', mv2) = normalise vars iv2
Sqrt iv
| Just v <- mv
-> let r = sqrt v
in ( ival r, Just r )
| otherwise
-> ( expr', Nothing )
where
expr' = Sqrt iv'
(iv', mv) = normalise vars iv
Cos iv
| Just v <- mv
-> let r = cos v
@ -298,7 +349,7 @@ normalise expr = case expr of
-> ( expr', Nothing )
where
expr' = Cos iv'
(iv', mv) = normalise iv
(iv', mv) = normalise vars iv
Sin iv
| Just v <- mv
-> let r = sin v
@ -307,7 +358,16 @@ normalise expr = case expr of
-> ( expr', Nothing )
where
expr' = Sin iv'
(iv', mv) = normalise iv
(iv', mv) = normalise vars iv
Atan iv
| Just v <- mv
-> let r = atan v
in ( ival r, Just r )
| otherwise
-> ( expr', Nothing )
where
expr' = Atan iv'
(iv', mv) = normalise vars iv
Pi -> ( Pi, Just pi )
Coerce1 iv
| Just v <- mv
@ -317,7 +377,7 @@ normalise expr = case expr of
-> ( expr', Nothing )
where
expr' = Coerce1 iv'
(iv', mv) = normalise iv
(iv', mv) = normalise vars iv
Coerce2 iv
| Just v <- mv
-> let r = coerce v
@ -326,7 +386,32 @@ normalise expr = case expr of
-> ( expr', Nothing )
where
expr' = Coerce2 iv'
(iv', mv) = normalise iv
(iv', mv) = normalise vars iv
EvalBez2 @v (Quadratic.Bezier p0 p1 p2) t
-> ( EvalBez2 (Quadratic.Bezier p0' p1' p2') t'
, Nothing
-- Quadratic.bezier
-- <$> (Quadratic.Bezier <$> m0 <*> m1 <*> m2)
-- <*> mt
)
where
(p0', m0) = normalise vars p0
(p1', m1) = normalise vars p1
(p2', m2) = normalise vars p2
(t' , mt) = normalise vars t
EvalBez3 @v (Cubic.Bezier p0 p1 p2 p3) t
-> ( EvalBez3 (Cubic.Bezier p0' p1' p2' p3') t'
, Nothing
-- Cubic.bezier
-- <$> (Cubic.Bezier <$> m0 <*> m1 <*> m2 <*> m3)
-- <*> mt
)
where
(p0', m0) = normalise vars p0
(p1', m1) = normalise vars p1
(p2', m2) = normalise vars p2
(p3', m3) = normalise vars p3
(t' , mt) = normalise vars t
Tabulate @n v
| Just xs <- sequence mxs
-> let r = tabulate ( xs ! )
@ -337,10 +422,10 @@ normalise expr = case expr of
expr' = Tabulate v'
v' :: Vec n ( AI Double )
mxs :: Vec n ( Maybe ( 𝕀 Double ) )
(v', mxs) = unzip $ fmap normalise v
(v', mxs) = unzip $ fmap ( normalise vars ) v
Index iv i
| Just r <- index_maybe iv' i
-> normalise r
-> normalise vars r
| Just v <- mv
-> let r = v `index` i
in ( ival r, Just r )
@ -348,13 +433,21 @@ normalise expr = case expr of
-> ( expr', Nothing )
where
expr' = Index iv' i
(iv', mv) = normalise iv
(iv', mv) = normalise vars iv
val :: Map String Double -> Value a -> a
val :: forall a. Typeable a => Map String Dynamic -> Value a -> a
val _ ( Val a ) = a
val vars ( Var v ) = vars Map.! v
val vars ( Var v ) =
case Map.lookup v vars of
Nothing -> error $ "Variable '" ++ v ++ "' not in environment."
Just d ->
case fromDynamic d of
Nothing -> error $ unlines [ "Incorrect type for variable '" ++ v ++ "'."
, " Expected: " ++ show ( typeRep $ Proxy @a )
, " Actual: " ++ show ( dynTypeRep d ) ]
Just a -> a
eval :: Map String Double -> AI a -> 𝕀 a
eval :: Typeable a => Map String Dynamic -> AI a -> 𝕀 a
eval vars = \ case
Pt v -> let f = val vars v in 𝕀 f f
Ival lo hi -> 𝕀 (val vars lo) (val vars hi)
@ -369,11 +462,15 @@ eval vars = \ case
Recip iv -> recip $ eval vars iv
Scale k v -> unT $ eval vars k *^ T ( eval vars v )
Cross iv1 iv2 -> T ( eval vars iv1 ) × T ( eval vars iv2 )
Sqrt iv -> sqrt $ eval vars iv
Cos iv -> cos $ eval vars iv
Sin iv -> sin $ eval vars iv
Atan iv -> atan $ eval vars iv
Pi -> pi
Coerce1 a -> coerce $ eval vars a
Coerce2 a -> coerce $ eval vars a
EvalBez2 bez t -> error "TODO" -- Quadratic.bezier (fmap (eval vars) bez) (eval vars t)
EvalBez3 bez t -> error "TODO" -- Cubic.bezier (fmap (eval vars) bez) (eval vars t)
Tabulate v -> let !w = fmap ( eval vars ) v
in tabulate ( w ! )
Index v i -> let !w = eval vars v
@ -390,7 +487,7 @@ instance Prelude.Num ( AI Double ) where
instance Prelude.Fractional ( AI Double ) where
recip = Recip
(/) = error "No division for abstract intervals"
a / b = a :*: Recip b
fromRational = Pt . Val . fromRational
instance Ring ( AI Double ) where
@ -404,10 +501,14 @@ deriving via ViaPrelude ( AI Double )
deriving via ViaPrelude ( AI Double )
instance Field ( AI Double )
instance Floating ( AI Double ) where
sqrt = Sqrt
instance Transcendental ( AI Double ) where
pi = Pi
cos = Cos
sin = Sin
atan = Atan
class ( Eq v, Show v, Module ( 𝕀 Double ) ( T ( 𝕀 v ) ), Group ( T ( 𝕀 v ) ) )
=> IsVec v

View file

@ -1,6 +1,7 @@
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE RebindableSyntax #-}
{-# LANGUAGE ScopedTypeVariables #-}
module MetaBrush.Asset.Brushes
@ -12,9 +13,9 @@ module MetaBrush.Asset.Brushes
-- base
import Prelude
hiding
( Num(..), Floating(..) )
( Num(..), Floating(..), (^), (/), fromInteger, fromRational )
import GHC.Exts
( Proxy# )
( Proxy#, fromString )
-- containers
import qualified Data.Sequence as Seq
@ -117,16 +118,16 @@ circleBrush _ mkI =
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 )
= ( x * r ) *^ e_x
^+^ ( y * r ) *^ e_y
mkPt x y
= ( r `scaledBy` x ) *^ e_x
^+^ ( r `scaledBy` y ) *^ e_y
in circleSpline @i @k @( Record CircleBrushFields ) @( 2 ) mkPt
where
e_x, e_y :: D k irec ( I i ( 2 ) )
e_x = pure $ mkI $ 2 1 0
e_y = pure $ mkI $ 2 0 1
kon = konst @( I i Double ) @k @irec . mkI
scaledBy d x = fmap ( mkI x * ) d
{-# INLINEABLE circleBrush #-}
ellipseBrush :: forall {t} (i :: t) k irec
@ -152,14 +153,37 @@ ellipseBrush _ mkI =
b = runD ( var @_ @k ( Fin 2 ) ) params
phi = runD ( var @_ @k ( Fin 3 ) ) params
mkPt :: Double -> Double -> D k irec ( I i ( 2 ) )
mkPt ( kon -> x ) ( kon -> y )
= ( x * a * cos phi - y * b * sin phi ) *^ e_x
^+^ ( y * b * cos phi + x * a * sin phi ) *^ e_y
mkPt x y
= let !x' = a `scaledBy` x
!y' = b `scaledBy` y
-- {-
in
( x' * cos phi - y' * sin phi ) *^ e_x
^+^ ( y' * cos phi + x' * sin phi ) *^ e_y
-- -}
{-
r = sqrt ( x' ^ 2 + y' ^ 2 )
arctgt = atan ( y' / x' )
-- a and b are always strictly positive, so we can compute
-- the quadrant using only x and y, which are constants.
!theta
| x > 0
= arctgt
| x < 0
= if y >= 0 then arctgt + pi else arctgt - pi
| otherwise
= if y >= 0 then 0.5 * pi else -0.5 * pi
!phi' = phi + theta
in
( r * cos phi' ) *^ e_x
^+^ ( r * sin phi' ) *^ e_y
-}
in circleSpline @i @k @( Record EllipseBrushFields ) @( 2 ) mkPt
where
e_x, e_y :: D k irec ( I i ( 2 ) )
e_x = pure $ mkI $ 2 1 0
e_y = pure $ mkI $ 2 0 1
kon = konst @( I i Double ) @k @irec . mkI
scaledBy d x = fmap ( mkI x * ) d
{-# INLINEABLE ellipseBrush #-}

View file

@ -25,15 +25,14 @@ module Math.Algebra.Dual
) where
-- base
import Prelude hiding ( Num(..), Floating(..), (^) )
import Control.Applicative
( liftA2 )
import Prelude
hiding
( Num(..), Floating(..)
, (^), recip, fromInteger, fromRational )
import Data.Coerce
( coerce )
import Data.Kind
( Type )
import Data.Monoid
( Ap(..) )
import GHC.TypeNats
( Nat )
@ -238,24 +237,37 @@ deriving via ApAp2 3 ( 4 ) r
--------------------------------------------------------------------------------
-- Ring instances.
-- TODO: should also implement a power operation specifically,
-- as this is important for interval arithmethic.
d1pow :: Ring r => Word -> r -> D1𝔸1 r
d1pow i x =
let !j = fromInteger $ fromIntegral i
in D11
( x ^ i )
( T $ j * x ^ (i - 1) )
{-# INLINE d1pow #-}
d2pow :: Ring r => Word -> r -> D2𝔸1 r
d2pow i x =
let !j = fromInteger $ fromIntegral i
in D21
( x ^ i )
( T $ j * x ^ (i - 1) )
( T $ j * (j - 1) * x ^ ( i - 2 ) )
{-# INLINE d2pow #-}
d3pow :: Ring r => Word -> r -> D3𝔸1 r
d3pow i x =
let !j = fromInteger $ fromIntegral i
in D31
( x ^ i )
( T $ j * x ^ (i - 1) )
( T $ j * (j - 1) * x ^ ( i - 2 ) )
( T $ j * (j - 1) * (j - 2) * x ^ ( i - 3 ) )
{-# INLINE d3pow #-}
deriving newtype instance Ring r => Ring ( D𝔸0 r )
--instance Ring r => Ring ( D1𝔸1 r ) where
-- !dr1 * !dr2 =
-- let
-- o :: r
-- o = fromInteger 0
-- p :: r -> r -> r
-- p = (+) @r
-- m :: r -> r -> r
-- m = (*) @r
-- in
-- $$( prodRuleQ
-- [|| o ||] [|| p ||] [|| m ||]
-- [|| dr1 ||] [|| dr2 ||] )
--instance Ring r => Ring ( D1𝔸2 r ) where
--instance Ring r => Ring ( D1𝔸3 r ) where
--instance Ring r => Ring ( D1𝔸4 r ) where
instance Ring r => Ring ( D2𝔸1 r ) where
!dr1 * !dr2 =
@ -270,39 +282,18 @@ instance Ring r => Ring ( D2𝔸1 r ) where
$$( prodRuleQ
[|| o ||] [|| p ||] [|| m ||]
[|| dr1 ||] [|| dr2 ||] )
df ^ i =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d2pow @r i ( _D21_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Ring ( D2𝔸1 Double ) #-}
{-# SPECIALISE instance Ring ( D2𝔸1 ( 𝕀 Double ) ) #-}
instance Ring r => Ring ( D3𝔸1 r ) where
!dr1 * !dr2 =
let
o :: r
!o = fromInteger 0
p :: r -> r -> r
!p = (+) @r
m :: r -> r -> r
!m = (*) @r
in
$$( prodRuleQ
[|| o ||] [|| p ||] [|| m ||]
[|| dr1 ||] [|| dr2 ||] )
{-# SPECIALISE instance Ring ( D3𝔸1 Double ) #-}
{-# SPECIALISE instance Ring ( D3𝔸1 ( 𝕀 Double ) ) #-}
--instance Ring r => Ring ( D1𝔸2 r ) where
-- !dr1 * !dr2 =
-- let
-- o :: r
-- o = fromInteger 0
-- p :: r -> r -> r
-- p = (+) @r
-- m :: r -> r -> r
-- m = (*) @r
-- in
-- $$( prodRuleQ
-- [|| o ||] [|| p ||] [|| m ||]
-- [|| dr1 ||] [|| dr2 ||] )
instance Ring r => Ring ( D2𝔸2 r ) where
!dr1 * !dr2 =
let
@ -316,39 +307,18 @@ instance Ring r => Ring ( D2𝔸2 r ) where
$$( prodRuleQ
[|| o ||] [|| p ||] [|| m ||]
[|| dr1 ||] [|| dr2 ||] )
df ^ i =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d2pow @r i ( _D22_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Ring ( D2𝔸2 Double ) #-}
{-# SPECIALISE instance Ring ( D2𝔸2 ( 𝕀 Double ) ) #-}
instance Ring r => Ring ( D3𝔸2 r ) where
!dr1 * !dr2 =
let
o :: r
!o = fromInteger 0
p :: r -> r -> r
!p = (+) @r
m :: r -> r -> r
!m = (*) @r
in
$$( prodRuleQ
[|| o ||] [|| p ||] [|| m ||]
[|| dr1 ||] [|| dr2 ||] )
{-# SPECIALISE instance Ring ( D3𝔸2 Double ) #-}
{-# SPECIALISE instance Ring ( D3𝔸2 ( 𝕀 Double ) ) #-}
--instance Ring r => Ring ( D1𝔸3 r ) where
-- !dr1 * !dr2 =
-- let
-- o :: r
-- o = fromInteger 0
-- p :: r -> r -> r
-- p = (+) @r
-- m :: r -> r -> r
-- m = (*) @r
-- in
-- $$( prodRuleQ
-- [|| o ||] [|| p ||] [|| m ||]
-- [|| dr1 ||] [|| dr2 ||] )
instance Ring r => Ring ( D2𝔸3 r ) where
!dr1 * !dr2 =
let
@ -362,39 +332,18 @@ instance Ring r => Ring ( D2𝔸3 r ) where
$$( prodRuleQ
[|| o ||] [|| p ||] [|| m ||]
[|| dr1 ||] [|| dr2 ||] )
df ^ i =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d2pow @r i ( _D23_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Ring ( D2𝔸3 Double ) #-}
{-# SPECIALISE instance Ring ( D2𝔸3 ( 𝕀 Double ) ) #-}
instance Ring r => Ring ( D3𝔸3 r ) where
!dr1 * !dr2 =
let
o :: r
!o = fromInteger 0
p :: r -> r -> r
!p = (+) @r
m :: r -> r -> r
!m = (*) @r
in
$$( prodRuleQ
[|| o ||] [|| p ||] [|| m ||]
[|| dr1 ||] [|| dr2 ||] )
{-# SPECIALISE instance Ring ( D3𝔸3 Double ) #-}
{-# SPECIALISE instance Ring ( D3𝔸3 ( 𝕀 Double ) ) #-}
--instance Ring r => Ring ( D1𝔸4 r ) where
-- !dr1 * !dr2 =
-- let
-- o :: r
-- o = fromInteger 0
-- p :: r -> r -> r
-- p = (+) @r
-- m :: r -> r -> r
-- m = (*) @r
-- in
-- $$( prodRuleQ
-- [|| o ||] [|| p ||] [|| m ||]
-- [|| dr1 ||] [|| dr2 ||] )
instance Ring r => Ring ( D2𝔸4 r ) where
!dr1 * !dr2 =
let
@ -408,9 +357,93 @@ instance Ring r => Ring ( D2𝔸4 r ) where
$$( prodRuleQ
[|| o ||] [|| p ||] [|| m ||]
[|| dr1 ||] [|| dr2 ||] )
df ^ i =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d2pow @r i ( _D24_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Ring ( D2𝔸4 Double ) #-}
{-# SPECIALISE instance Ring ( D2𝔸4 ( 𝕀 Double ) ) #-}
instance Ring r => Ring ( D3𝔸1 r ) where
!dr1 * !dr2 =
let
o :: r
!o = fromInteger 0
p :: r -> r -> r
!p = (+) @r
m :: r -> r -> r
!m = (*) @r
in
$$( prodRuleQ
[|| o ||] [|| p ||] [|| m ||]
[|| dr1 ||] [|| dr2 ||] )
df ^ i =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d3pow @r i ( _D31_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Ring ( D3𝔸1 Double ) #-}
{-# SPECIALISE instance Ring ( D3𝔸1 ( 𝕀 Double ) ) #-}
instance Ring r => Ring ( D3𝔸2 r ) where
!dr1 * !dr2 =
let
o :: r
!o = fromInteger 0
p :: r -> r -> r
!p = (+) @r
m :: r -> r -> r
!m = (*) @r
in
$$( prodRuleQ
[|| o ||] [|| p ||] [|| m ||]
[|| dr1 ||] [|| dr2 ||] )
df ^ i =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d3pow @r i ( _D32_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Ring ( D3𝔸2 Double ) #-}
{-# SPECIALISE instance Ring ( D3𝔸2 ( 𝕀 Double ) ) #-}
instance Ring r => Ring ( D3𝔸3 r ) where
!dr1 * !dr2 =
let
o :: r
!o = fromInteger 0
p :: r -> r -> r
!p = (+) @r
m :: r -> r -> r
!m = (*) @r
in
$$( prodRuleQ
[|| o ||] [|| p ||] [|| m ||]
[|| dr1 ||] [|| dr2 ||] )
df ^ i =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d3pow @r i ( _D33_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Ring ( D3𝔸3 Double ) #-}
{-# SPECIALISE instance Ring ( D3𝔸3 ( 𝕀 Double ) ) #-}
instance Ring r => Ring ( D3𝔸4 r ) where
!dr1 * !dr2 =
let
@ -424,11 +457,36 @@ instance Ring r => Ring ( D3𝔸4 r ) where
$$( prodRuleQ
[|| o ||] [|| p ||] [|| m ||]
[|| dr1 ||] [|| dr2 ||] )
df ^ i =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d2pow @r i ( _D34_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Ring ( D3𝔸4 Double ) #-}
{-# SPECIALISE instance Ring ( D3𝔸4 ( 𝕀 Double ) ) #-}
--------------------------------------------------------------------------------
-- Field & transcendental instances
-- Field, floating & transcendental instances
d1recip :: Field r => r -> D1𝔸1 r
d1recip x =
let !d = recip x
in D11 d ( T $ -1 * d ^ 2 )
{-# INLINE d1recip #-}
d2recip :: Field r => r -> D2𝔸1 r
d2recip x =
let !d = recip x
in D21 d ( T $ -1 * d ^ 2 ) ( T $ 2 * d ^ 3 )
{-# INLINE d2recip #-}
d3recip :: Field r => r -> D3𝔸1 r
d3recip x =
let !d = recip x
in D31 d ( T $ -1 * d ^ 2 ) ( T $ 2 * d ^ 3 ) ( T $ -6 * d ^ 4 )
{-# INLINE d3recip #-}
deriving newtype instance Field r => Field ( D𝔸0 r )
@ -437,16 +495,232 @@ deriving newtype instance Field r => Field ( D𝔸0 r )
--instance Field r => Field ( D1𝔸3 r ) where
--instance Field r => Field ( D1𝔸4 r ) where
instance Field r => Field ( D2𝔸1 r ) where
fromRational q = konst @Double @2 @( 1 ) ( fromRational q )
recip df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d2recip @r ( _D21_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Field ( D2𝔸1 Double ) #-}
{-# SPECIALISE instance Field ( D2𝔸1 ( 𝕀 Double ) ) #-}
instance Field r => Field ( D2𝔸2 r ) where
fromRational q = konst @Double @2 @( 2 ) ( fromRational q )
recip df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d2recip @r ( _D22_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Field ( D2𝔸2 Double ) #-}
{-# SPECIALISE instance Field ( D2𝔸2 ( 𝕀 Double ) ) #-}
instance Field r => Field ( D2𝔸3 r ) where
fromRational q = konst @Double @2 @( 3 ) ( fromRational q )
recip df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d2recip @r ( _D23_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Field ( D2𝔸3 Double ) #-}
{-# SPECIALISE instance Field ( D2𝔸3 ( 𝕀 Double ) ) #-}
instance Field r => Field ( D2𝔸4 r ) where
fromRational q = konst @Double @2 @( 4 ) ( fromRational q )
recip df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d2recip @r ( _D24_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Field ( D2𝔸4 Double ) #-}
{-# SPECIALISE instance Field ( D2𝔸4 ( 𝕀 Double ) ) #-}
instance Field r => Field ( D3𝔸1 r ) where
fromRational q = konst @Double @3 @( 1 ) ( fromRational q )
recip df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d3recip @r ( _D31_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Field ( D3𝔸1 Double ) #-}
{-# SPECIALISE instance Field ( D3𝔸1 ( 𝕀 Double ) ) #-}
instance Field r => Field ( D3𝔸2 r ) where
fromRational q = konst @Double @3 @( 2 ) ( fromRational q )
recip df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d3recip @r ( _D32_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Field ( D3𝔸2 Double ) #-}
{-# SPECIALISE instance Field ( D3𝔸2 ( 𝕀 Double ) ) #-}
instance Field r => Field ( D3𝔸3 r ) where
fromRational q = konst @Double @3 @( 3 ) ( fromRational q )
recip df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d3recip @r ( _D33_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Field ( D3𝔸3 Double ) #-}
{-# SPECIALISE instance Field ( D3𝔸3 ( 𝕀 Double ) ) #-}
instance Field r => Field ( D3𝔸4 r ) where
-- TODO
fromRational q = konst @Double @3 @( 4 ) ( fromRational q )
recip df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d3recip @r ( _D34_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Field ( D3𝔸4 Double ) #-}
{-# SPECIALISE instance Field ( D3𝔸4 ( 𝕀 Double ) ) #-}
d1sin, d1cos :: Transcendental r => r -> D1𝔸1 r
d1sqrt :: Floating r => r -> D1𝔸1 r
d1sqrt x =
let !v = sqrt x
!d = recip v
in D11 v ( T $ 0.5 * d )
{-# INLINE d1sqrt #-}
d2sqrt :: Floating r => r -> D2𝔸1 r
d2sqrt x =
let !v = sqrt x
!d = recip v
in D21 v ( T $ 0.5 * d ) ( T $ -0.25 * d ^ 2 )
{-# INLINE d2sqrt #-}
d3sqrt :: Floating r => r -> D3𝔸1 r
d3sqrt x =
let !v = sqrt x
!d = recip v
in D31 v ( T $ 0.5 * d ) ( T $ -0.25 * d ^ 2 ) ( T $ 0.375 * d ^ 3 )
{-# INLINE d3sqrt #-}
deriving newtype instance Floating r => Floating ( D𝔸0 r )
--instance Floating r => Floating ( D1𝔸1 r ) where
--instance Floating r => Floating ( D1𝔸2 r ) where
--instance Floating r => Floating ( D1𝔸3 r ) where
--instance Floating r => Floating ( D1𝔸4 r ) where
instance Floating r => Floating ( D2𝔸1 r ) where
sqrt df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d2sqrt @r ( _D21_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Floating ( D2𝔸1 Double ) #-}
{-# SPECIALISE instance Floating ( D2𝔸1 ( 𝕀 Double ) ) #-}
instance Floating r => Floating ( D2𝔸2 r ) where
sqrt df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d2sqrt @r ( _D22_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Floating ( D2𝔸2 Double ) #-}
{-# SPECIALISE instance Floating ( D2𝔸2 ( 𝕀 Double ) ) #-}
instance Floating r => Floating ( D2𝔸3 r ) where
sqrt df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d2sqrt @r ( _D23_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Floating ( D2𝔸3 Double ) #-}
{-# SPECIALISE instance Floating ( D2𝔸3 ( 𝕀 Double ) ) #-}
instance Floating r => Floating ( D2𝔸4 r ) where
sqrt df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d2sqrt @r ( _D24_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Floating ( D2𝔸4 Double ) #-}
{-# SPECIALISE instance Floating ( D2𝔸4 ( 𝕀 Double ) ) #-}
instance Floating r => Floating ( D3𝔸1 r ) where
sqrt df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d3sqrt @r ( _D31_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Floating ( D3𝔸1 Double ) #-}
{-# SPECIALISE instance Floating ( D3𝔸1 ( 𝕀 Double ) ) #-}
instance Floating r => Floating ( D3𝔸2 r ) where
sqrt df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d3sqrt @r ( _D32_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Floating ( D3𝔸2 Double ) #-}
{-# SPECIALISE instance Floating ( D3𝔸2 ( 𝕀 Double ) ) #-}
instance Floating r => Floating ( D3𝔸3 r ) where
sqrt df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d3sqrt @r ( _D33_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Floating ( D3𝔸3 Double ) #-}
{-# SPECIALISE instance Floating ( D3𝔸3 ( 𝕀 Double ) ) #-}
instance Floating r => Floating ( D3𝔸4 r ) where
sqrt df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d3sqrt @r ( _D34_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Floating ( D3𝔸4 Double ) #-}
{-# SPECIALISE instance Floating ( D3𝔸4 ( 𝕀 Double ) ) #-}
d1sin, d1cos, d1atan :: Transcendental r => r -> D1𝔸1 r
d1sin x =
let !s = sin x
!c = cos x
@ -455,8 +729,16 @@ d1cos x =
let !s = sin x
!c = cos x
in D11 c ( T -s )
d1atan x
= let !d = recip $ 1 + x ^ 2
in D11
( atan x )
( T $ d )
{-# INLINE d1sin #-}
{-# INLINE d1cos #-}
{-# INLINE d1atan #-}
d2sin, d2cos :: Transcendental r => r -> D2𝔸1 r
d2sin, d2cos, d2atan :: Transcendental r => r -> D2𝔸1 r
d2sin x =
let !s = sin x
!c = cos x
@ -465,8 +747,17 @@ d2cos x =
let !s = sin x
!c = cos x
in D21 c ( T -s ) ( T -c )
d2atan x
= let !d = recip $ 1 + x ^ 2
in D21
( atan x )
( T $ d )
( T $ -2 * x * d ^ 2 )
{-# INLINE d2sin #-}
{-# INLINE d2cos #-}
{-# INLINE d2atan #-}
d3sin, d3cos :: Transcendental r => r -> D3𝔸1 r
d3sin, d3cos, d3atan :: Transcendental r => r -> D3𝔸1 r
d3sin x =
let !s = sin x
!c = cos x
@ -475,6 +766,16 @@ d3cos x =
let !s = sin x
!c = cos x
in D31 c ( T -s ) ( T -c ) ( T s )
d3atan x
= let !d = recip $ 1 + x ^ 2
in D31
( atan x )
( T $ d )
( T $ -2 * x * d ^ 2 )
( T $ ( -2 + 6 * x ^ 2 ) * d ^ 3 )
{-# INLINE d3sin #-}
{-# INLINE d3cos #-}
{-# INLINE d3atan #-}
deriving newtype instance Transcendental r => Transcendental ( D𝔸0 r )
--instance Transcendental r => Transcendental ( D1𝔸1 r ) where
@ -501,6 +802,17 @@ instance Transcendental r => Transcendental ( D2𝔸1 r ) where
dg = d2cos @r ( _D21_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
atan df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d2atan @r ( _D21_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Transcendental ( D2𝔸1 Double ) #-}
{-# SPECIALISE instance Transcendental ( D2𝔸1 ( 𝕀 Double ) ) #-}
instance Transcendental r => Transcendental ( D2𝔸2 r ) where
pi = konst @Double @2 @( 2 ) pi
sin df =
@ -521,6 +833,17 @@ instance Transcendental r => Transcendental ( D2𝔸2 r ) where
dg = d2cos @r ( _D22_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
atan df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d2atan @r ( _D22_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Transcendental ( D2𝔸2 Double ) #-}
{-# SPECIALISE instance Transcendental ( D2𝔸2 ( 𝕀 Double ) ) #-}
instance Transcendental r => Transcendental ( D2𝔸3 r ) where
pi = konst @Double @2 @( 3 ) pi
sin df =
@ -541,6 +864,17 @@ instance Transcendental r => Transcendental ( D2𝔸3 r ) where
dg = d2cos @r ( _D23_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
atan df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d2atan @r ( _D23_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Transcendental ( D2𝔸3 Double ) #-}
{-# SPECIALISE instance Transcendental ( D2𝔸3 ( 𝕀 Double ) ) #-}
instance Transcendental r => Transcendental ( D2𝔸4 r ) where
pi = konst @Double @2 @( 4 ) pi
@ -562,6 +896,17 @@ instance Transcendental r => Transcendental ( D2𝔸4 r ) where
dg = d2cos @r ( _D24_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
atan df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d2atan @r ( _D24_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Transcendental ( D2𝔸4 Double ) #-}
{-# SPECIALISE instance Transcendental ( D2𝔸4 ( 𝕀 Double ) ) #-}
instance Transcendental r => Transcendental ( D3𝔸1 r ) where
pi = konst @Double @3 @( 1 ) pi
@ -583,6 +928,17 @@ instance Transcendental r => Transcendental ( D3𝔸1 r ) where
dg = d3cos @r ( _D31_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
atan df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d3atan @r ( _D31_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Transcendental ( D3𝔸1 Double ) #-}
{-# SPECIALISE instance Transcendental ( D3𝔸1 ( 𝕀 Double ) ) #-}
instance Transcendental r => Transcendental ( D3𝔸2 r ) where
pi = konst @Double @3 @( 2 ) pi
@ -604,6 +960,17 @@ instance Transcendental r => Transcendental ( D3𝔸2 r ) where
dg = d3cos @r ( _D32_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
atan df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d3atan @r ( _D32_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Transcendental ( D3𝔸2 Double ) #-}
{-# SPECIALISE instance Transcendental ( D3𝔸2 ( 𝕀 Double ) ) #-}
instance Transcendental r => Transcendental ( D3𝔸3 r ) where
pi = konst @Double @3 @( 3 ) pi
@ -625,6 +992,17 @@ instance Transcendental r => Transcendental ( D3𝔸3 r ) where
dg = d3cos @r ( _D33_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
atan df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d3atan @r ( _D33_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Transcendental ( D3𝔸3 Double ) #-}
{-# SPECIALISE instance Transcendental ( D3𝔸3 ( 𝕀 Double ) ) #-}
instance Transcendental r => Transcendental ( D3𝔸4 r ) where
pi = konst @Double @3 @( 4 ) pi
@ -646,6 +1024,17 @@ instance Transcendental r => Transcendental ( D3𝔸4 r ) where
dg = d3cos @r ( _D34_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
atan df =
let
fromInt = fromInteger @r
add = (+) @r
times = (*) @r
pow = (^) @r
dg = d3atan @r ( _D34_v df )
in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||]
[|| df ||] [|| dg ||] )
{-# SPECIALISE instance Transcendental ( D3𝔸4 Double ) #-}
{-# SPECIALISE instance Transcendental ( D3𝔸4 ( 𝕀 Double ) ) #-}
--------------------------------------------------------------------------------
-- HasChainRule instances.

View file

@ -22,8 +22,6 @@ module Math.Bezier.Stroke
-- base
import Control.Arrow
( first, (***) )
import Control.Applicative
( Applicative(..) )
import Control.Monad
( unless )
import Control.Monad.ST
@ -31,7 +29,7 @@ import Control.Monad.ST
import Data.Bifunctor
( Bifunctor(bimap) )
import Data.Coerce
( Coercible, coerce )
( coerce )
import Data.Foldable
( for_ )
import Data.Functor.Identity
@ -1273,7 +1271,7 @@ intervalNewtonGS precondMethod minWidth eqs =
, 𝛿E𝛿sdcdt = D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) }
<- ( eqs t `Seq.index` i ) s
, StrokeDatum { 𝛿E𝛿sdcdt = D12 ( T f_mid ) ( T ( T f_t_mid ) ) ( T ( T f_s_mid ) ) }
, 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 | inf ee < 1 0
, sup ee > 1 0

View file

@ -19,8 +19,6 @@ module Math.Interval
-- base
import Prelude hiding ( Num(..) )
import Data.List.NonEmpty
( NonEmpty(..) )
-- acts
import Data.Act

View file

@ -116,12 +116,24 @@ instance Prelude.Fractional ( 𝕀 Double ) where
| otherwise
= error "BAD interval recip; should use extendedRecip instead"
-- #endif
_ / _ = error "TODO: interval division is not implemented"
𝕀 x_lo x_hi / 𝕀 y_lo y_hi
-- #ifdef ASSERTS
| y_lo >= 0 || y_hi <= 0
-- #endif
= 𝕀 ( fst $ divI x_lo y_hi ) ( snd $ divI x_hi y_lo )
-- #ifdef ASSERTS
| otherwise
= error "BAD interval division; should use extendedRecip instead"
-- #endif
instance Floating ( 𝕀 Double ) where
sqrt = withHW Prelude.sqrt
instance Transcendental ( 𝕀 Double ) where
pi = 𝕀 3.141592653589793 3.1415926535897936
cos = withHW Prelude.cos
sin = withHW Prelude.sin
atan = withHW Prelude.atan
{-# INLINE withHW #-}
withHW :: (Interval.Interval a -> Interval.Interval b) -> 𝕀 a -> 𝕀 b

View file

@ -96,7 +96,6 @@ instance {-# INCOHERENT #-} Show v => Show ( T v ) where
showsPrec p ( T v )
= showParen ( p >= 11 )
$ showString "T"
. showSpace
. showsPrec 0 v
instance Applicative T where

View file

@ -19,8 +19,6 @@ module Math.Module
where
-- base
import Control.Applicative
( liftA2 )
import Control.Monad
( guard )
import Data.Coerce

View file

@ -14,8 +14,6 @@ module Math.Ring
import Prelude ( Num, Fractional )
import Prelude hiding ( Num(..), Fractional(..), Floating(..) )
import qualified Prelude
import Control.Applicative
( liftA2 )
import Data.Coerce
( coerce )
import Data.Monoid
@ -64,10 +62,11 @@ class Ring a => Field a where
class Field a => Floating a where
sqrt :: a -> a
class Field a => Transcendental a where
class Floating a => Transcendental a where
pi :: a
cos :: a -> a
sin :: a -> a
atan :: a -> a
--------------------------------------------------------------------------------
@ -136,6 +135,7 @@ instance Prelude.Floating a => Transcendental ( ViaPrelude a ) where
pi = coerce $ Prelude.pi @a
sin = coerce $ Prelude.sin @a
cos = coerce $ Prelude.cos @a
atan = coerce $ Prelude.atan @a
--------------------------------------------------------------------------------