diff --git a/src/metabrushes/MetaBrush/Assert.hs b/src/metabrushes/MetaBrush/Assert.hs index 1a05d1e..a9e4153 100644 --- a/src/metabrushes/MetaBrush/Assert.hs +++ b/src/metabrushes/MetaBrush/Assert.hs @@ -13,9 +13,9 @@ import Control.Exception -------------------------------------------------------------------------------- {-# INLINE assert #-} -assert :: String -> a -> a +assert :: Bool -> String -> a -> a #ifdef ASSERTS -assert message _ = throw ( AssertionFailed message ) +assert False message _ = throw ( AssertionFailed message ) #else -assert _ a = a +assert _ _ a = a #endif diff --git a/src/metabrushes/MetaBrush/Document/Draw.hs b/src/metabrushes/MetaBrush/Document/Draw.hs index 4ab9830..6ae752d 100644 --- a/src/metabrushes/MetaBrush/Document/Draw.hs +++ b/src/metabrushes/MetaBrush/Document/Draw.hs @@ -220,7 +220,7 @@ addToAnchor anchor newSpline = over ( field' @"documentContent" . field' @"strok setBrushData = set ( field @"brushParams" ) ( brushParams ( splineEnd prevSpline ) ) in prevSpline <> fmap setBrushData newSpline | otherwise - = assert ( "addToAnchor: trying to add to closed spline " <> show strokeUnique ) + = assert False ( "addToAnchor: trying to add to closed spline " <> show strokeUnique ) prevSpline -- should never add to a closed spline = overStrokeSpline updateSpline stroke | otherwise diff --git a/src/splines/Math/Bezier/Cubic.hs b/src/splines/Math/Bezier/Cubic.hs index 6bf3c50..e0f8a58 100644 --- a/src/splines/Math/Bezier/Cubic.hs +++ b/src/splines/Math/Bezier/Cubic.hs @@ -6,10 +6,12 @@ module Math.Bezier.Cubic ( Bezier(..) , fromQuadratic , bezier, bezier', bezier'', bezier''' + , derivative , curvature, squaredCurvature, signedCurvature , subdivide , ddist, closestPoint , drag, selfIntersectionParameters + , extrema ) where @@ -105,6 +107,11 @@ bezier ( Bezier {..} ) t = ( Quadratic.bezier @v ( Quadratic.Bezier p0 p1 p2 ) t ) ( Quadratic.bezier @v ( Quadratic.Bezier p1 p2 p3 ) t ) +-- | The derivative of a Cubic Bézier curve, as a quadratic Bézier curve. +derivative :: ( Group v, Module r v ) => Bezier v -> Quadratic.Bezier v +derivative ( Bezier {..} ) = ( Ring.fromInteger 3 *^ ) + <$> Quadratic.Bezier ( p0 --> p1 ) ( p1 --> p2 ) ( p2 --> p3 ) + -- | Derivative of a cubic Bézier curve. bezier' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v bezier' ( Bezier {..} ) @@ -255,3 +262,11 @@ selfIntersectionParameters ( Bezier {..} ) = solveQuadratic c0 c1 c2 c0 = f2 c1 = f3 - 2 * f2 c2 = f1 + f2 - f3 + +-- | Extremal values of the Bézier parameter for a cubic Bézier curve. +extrema :: RealFloat r => Bezier r -> [ r ] +extrema ( Bezier {..} ) = solveQuadratic c b a + where + a = p3 - 3 * p2 + 3 * p1 - p0 + b = 2 * ( p0 - 2 * p1 + p2 ) + c = p1 - p0 diff --git a/src/splines/Math/Bezier/Quadratic.hs b/src/splines/Math/Bezier/Quadratic.hs index 3da728d..f56493c 100644 --- a/src/splines/Math/Bezier/Quadratic.hs +++ b/src/splines/Math/Bezier/Quadratic.hs @@ -9,6 +9,7 @@ module Math.Bezier.Quadratic , subdivide , ddist, closestPoint , interpolate + , extrema ) where @@ -188,3 +189,9 @@ interpolate p0 p2 t q = Bezier {..} p1 = ( ( 0.5 * ( t - 1 ) / t ) *^ ( q --> p0 :: v ) ^+^ ( 0.5 * t / ( t - 1 ) ) *^ ( q --> p2 :: v ) ) • q + +-- | Extremal values of the Bézier parameter for a quadratic Bézier curve. +extrema :: Fractional r => Bezier r -> [ r ] +extrema ( Bezier {..} ) = [ t ] + where + t = ( p0 - p1 ) / ( p0 - 2 * p1 + p2 ) diff --git a/src/splines/Math/Bezier/Stroke.hs b/src/splines/Math/Bezier/Stroke.hs index 87750ed..27de833 100644 --- a/src/splines/Math/Bezier/Stroke.hs +++ b/src/splines/Math/Bezier/Stroke.hs @@ -235,6 +235,7 @@ computeStrokeOutline :: , HasChainRule Double 2 brushParams , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 brushParams ) , Traversable ( D 2 brushParams ) + , Representable Double usedParams -- Debugging. , Show ptData, Show brushParams @@ -513,6 +514,9 @@ outlineFunction , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 brushParams ) , Traversable ( D 2 brushParams ) + -- Computing AABBs + , Representable Double usedParams + -- Debugging. , Show ptData, Show brushParams ) @@ -600,13 +604,16 @@ outlineFunction ptParams toBrushParams brushFromParams = \ sp0 crv -> pathAndUsedParams :: forall k i arr crvData ptData usedParams . ( HasType ( ℝ 2 ) ptData - , CurveOrder k + , HasBézier k i , arr ~ C 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 ) ) , Module ( I i Double ) ( T ( I i usedParams ) ) , Torsor ( T ( I i usedParams ) ) ( I i usedParams ) + , Module Double ( T usedParams ) + , Representable Double usedParams + , Torsor ( T usedParams ) usedParams ) => ( I i ( ℝ 1 ) -> I i Double ) -> ( forall a. a -> I i a ) @@ -618,16 +625,16 @@ pathAndUsedParams co toI ptParams sp0 crv = case crv of LineTo { curveEnd = NextPoint sp1 } | let seg = Segment sp0 sp1 - -> ( line @k @i co ( fmap ( toI . coords ) seg ) - , line @k @i co ( fmap ( toI . ptParams ) seg ) ) + -> ( line @k @i @( ℝ 2 ) co ( fmap ( toI . coords ) seg ) + , line @k @i @usedParams co ( fmap ( toI . ptParams ) seg ) ) Bezier2To { controlPoint = sp1, curveEnd = NextPoint sp2 } | let bez2 = Quadratic.Bezier sp0 sp1 sp2 - -> ( bezier2 @k @i co ( fmap ( toI . coords ) bez2 ) - , bezier2 @k @i co ( fmap ( toI . ptParams ) bez2 ) ) + -> ( bezier2 @k @i @( ℝ 2 ) co ( fmap ( toI . coords ) bez2 ) + , bezier2 @k @i @usedParams co ( fmap ( toI . ptParams ) bez2 ) ) Bezier3To { controlPoint1 = sp1, controlPoint2 = sp2, curveEnd = NextPoint sp3 } | let bez3 = Cubic.Bezier sp0 sp1 sp2 sp3 - -> ( bezier3 @k @i co ( fmap ( toI . coords ) bez3 ) - , bezier3 @k @i co ( fmap ( toI . ptParams ) bez3 ) ) + -> ( bezier3 @k @i @( ℝ 2 ) co ( fmap ( toI . coords ) bez3 ) + , bezier3 @k @i @usedParams co ( fmap ( toI . ptParams ) bez3 ) ) ----------------------------------- -- Various utility functions @@ -923,7 +930,7 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) } splineCurveFns :: forall k i - . ( CurveOrder k + . ( HasBézier k i , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) , Module ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) ) ) @@ -940,14 +947,14 @@ splineCurveFns co spls where curveFn :: I i ( ℝ 2 ) -> Curve Open () ( I i ( ℝ 2 ) ) - -> ( C k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) ) + -> C k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) curveFn p0 = \case LineTo { curveEnd = NextPoint p1 } - -> line @k @i co $ Segment p0 p1 + -> line @k @i @( ℝ 2 ) co $ Segment p0 p1 Bezier2To { controlPoint = p1, curveEnd = NextPoint p2 } - -> bezier2 @k @i co $ Quadratic.Bezier p0 p1 p2 + -> bezier2 @k @i @( ℝ 2 ) co $ Quadratic.Bezier p0 p1 p2 Bezier3To { controlPoint1 = p1, controlPoint2 = p2, curveEnd = NextPoint p3 } - -> bezier3 @k @i co $ Cubic.Bezier p0 p1 p2 p3 + -> bezier3 @k @i @( ℝ 2 ) co $ Cubic.Bezier p0 p1 p2 p3 -- | Solve the envelope equations at a given point \( t = t_0 \), to find -- \( s_0 \) such that \( c(t_0, s_0) \) is on the envelope of the brush stroke. @@ -1072,7 +1079,8 @@ instance Applicative ZipSeq where {-# INLINE liftA2 #-} brushStrokeData :: forall k brushParams i arr - . ( CurveOrder k, arr ~ C k + . ( arr ~ C k + , HasBézier k i, HasEnvelopeEquation k , Differentiable k i brushParams , HasChainRule ( I i Double ) k ( I i ( ℝ 1 ) ) , Applicative ( D k ( ℝ 1 ) ) @@ -1084,7 +1092,6 @@ brushStrokeData :: forall k brushParams i arr , D k ( I i ( ℝ 2 ) ) ~ D k ( ℝ 2 ) , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) ) - -- , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) , Show brushParams ) => ( I i Double -> I i ( ℝ 1 ) ) diff --git a/src/splines/Math/Bezier/Stroke/EnvelopeEquation.hs b/src/splines/Math/Bezier/Stroke/EnvelopeEquation.hs index 594eb9c..5ce25e6 100644 --- a/src/splines/Math/Bezier/Stroke/EnvelopeEquation.hs +++ b/src/splines/Math/Bezier/Stroke/EnvelopeEquation.hs @@ -5,17 +5,16 @@ module Math.Bezier.Stroke.EnvelopeEquation ( StrokeDatum(..) - , CurveOrder(..) + , HasBézier(..) + , HasEnvelopeEquation(..) ) where -- base import Prelude hiding ( Num(..), (^) ) -import Data.Coerce - ( Coercible, coerce ) -import Data.Functor.Identity - ( Identity(..) ) import Data.Kind ( Type, Constraint ) +import Data.List.NonEmpty + ( NonEmpty(..) ) import GHC.TypeNats ( Nat, type (-) ) @@ -72,32 +71,40 @@ deriving stock instance Show ( StrokeDatum 3 𝕀 ) -------------------------------------------------------------------------------- -type CurveOrder :: Nat -> Constraint -class CurveOrder k where +type HasBézier :: forall {t}. Nat -> t -> Constraint +class HasBézier k i where -- | Linear interpolation, as a differentiable function. - line :: forall i b - . ( Module ( I i Double ) ( T b ), Torsor ( T b ) b + line :: forall b + . ( Module Double ( T b ), Torsor ( T b ) b + , Module ( I i Double ) ( T (I i b ) ), Torsor ( T ( I i b ) ) ( I i b ) , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) ) => ( I i ( ℝ 1 ) -> I i Double ) - -> Segment b -> C k ( I i ( ℝ 1 ) ) b + -> Segment ( I i b ) -> C k ( I i ( ℝ 1 ) ) ( I i b ) -- | A quadratic Bézier curve, as a differentiable function. - bezier2 :: forall i b - . ( Module ( I i Double ) ( T b ), Torsor ( T b ) b + bezier2 :: forall b + . ( Module Double ( T b ), Torsor ( T b ) b + , Module ( I i Double ) ( T (I i b ) ), Torsor ( T ( I i b ) ) ( I i b ) + , Representable Double b , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) ) => ( I i ( ℝ 1 ) -> I i Double ) - -> Quadratic.Bezier b -> C k ( I i ( ℝ 1 ) ) b + -> Quadratic.Bezier ( I i b ) -> C k ( I i ( ℝ 1 ) ) ( I i b ) -- | A cubic Bézier curve, as a differentiable function. - bezier3 :: forall i b - . ( Module ( I i Double ) ( T b ), Torsor ( T b ) b + bezier3 :: forall b + . ( Module Double ( T b ), Torsor ( T b ) b + , Module ( I i Double ) ( T (I i b ) ), Torsor ( T ( I i b ) ) ( I i b ) + , Representable Double b , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) ) => ( I i ( ℝ 1 ) -> I i Double ) - -> Cubic.Bezier b -> C k ( I i ( ℝ 1 ) ) b + -> Cubic.Bezier ( I i b ) -> C k ( I i ( ℝ 1 ) ) ( I i b ) + +type HasEnvelopeEquation :: Nat -> Constraint +class HasEnvelopeEquation k where uncurryD :: D k a ~ D k ( ℝ 1 ) => D k ( ℝ 1 ) ( C k a b ) -> a -> D k ( ℝ 2 ) b @@ -137,7 +144,7 @@ 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 +instance HasBézier 2 () where line co ( Segment a b :: Segment b ) = D \ ( co -> t ) -> D21 ( lerp @( T b ) t a b ) @@ -156,6 +163,8 @@ instance CurveOrder 2 where ( Cubic.bezier' bez t ) ( Cubic.bezier'' bez t ) +instance HasEnvelopeEquation 2 where + uncurryD = uncurryD2 brushStroke ( D21 p p_t p_tt ) ( D22 b b_t b_s b_tt b_ts b_ss ) = @@ -185,7 +194,7 @@ instance CurveOrder 2 where -- -- ∂E/∂s dc/dt = ∂E/∂s ∂c/∂t - ∂E/∂t ∂c/∂s. -instance CurveOrder 3 where +instance HasBézier 3 () where line co ( Segment a b :: Segment b ) = D \ ( co -> t ) -> @@ -208,6 +217,8 @@ instance CurveOrder 3 where ( Cubic.bezier'' bez t ) ( Cubic.bezier''' bez ) +instance HasEnvelopeEquation 3 where + uncurryD = uncurryD3 brushStroke @@ -246,3 +257,79 @@ instance CurveOrder 3 where ( T $ co ee_t ) ( T $ co ee_s ) ( T $ co ee_tt) ( T $ co ee_ts ) ( T $ co ee_ss ) , D12 𝛿E𝛿sdcdt ( T 𝛿E𝛿sdcdt_t ) ( T 𝛿E𝛿sdcdt_s ) ) + +instance HasBézier 3 𝕀 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 ( aabb bez ( `evaluateQuadratic` t ) ) + ( Quadratic.bezier' bez t ) + ( Quadratic.bezier'' bez ) + origin + + bezier3 co ( bez :: Cubic.Bezier b ) = + D \ ( co -> t ) -> + D31 ( aabb bez ( `evaluateCubic` t ) ) + ( T $ aabb ( fmap unT $ Cubic.derivative ( fmap T bez ) ) ( `evaluateQuadratic` t ) ) + ( Cubic.bezier'' bez t ) + ( Cubic.bezier''' bez ) + +{- Note [Computing Béziers over intervals] +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Here's how we evaluate a Bézier curve when the coefficients are intervals. +As we are using axis-aligned interval arithmetic, we reduce to the 1D situation. + +Now, the formulas for Bézier curves (with value of the t parameter in [0,1]) +are convex combinations + + b_0(t) [p_0,q_0] + b_1(t) [p_1,q_1] + ... + b_n(t) [p_n,q_n] + + -- Here b_1, ..., b_n are Bernstein polynomials. + +This means that the minimum value attained by the Bézier curve as we vary +both the time parameter and the values of the points within their respective +intervals will occur when we take + + b_0(t) p_0 + b_1(t) p_1 + ... + b_n(t) p_n + +and the maximum when we take + + b_0(t) q_0 + b_1(t) q_1 + ... + b_n(t) q_n + +For each of these, we can compute a minimum/maximum by computing an axis-aligned +bounding box for the Bézier curve. This is done by computing the derivative +with respect to t and root-finding. +-} + +-- | Evaluate a cubic Bézier curve, when both the coefficients and the +-- parameter are intervals. +evaluateCubic :: Cubic.Bezier ( 𝕀 Double ) -> 𝕀 Double -> 𝕀 Double +evaluateCubic bez t = + -- assert (inf t >= 0 && sup t <= 1) "evaluateCubic: t ⊊ [0,1]" $ -- Requires t ⊆ [0,1] + let inf_bez = fmap inf bez + sup_bez = fmap sup bez + mins = fmap (Cubic.bezier @( T Double ) inf_bez) + $ inf t :| ( sup t : filter ( inside t ) ( Cubic.extrema inf_bez ) ) + maxs = fmap (Cubic.bezier @( T Double ) sup_bez) + $ inf t :| ( sup t : filter ( inside t ) ( Cubic.extrema sup_bez ) ) + in 𝕀 ( minimum mins ) ( maximum maxs ) + +-- | Evaluate a quadratic Bézier curve, when both the coefficients and the +-- parameter are intervals. +evaluateQuadratic :: Quadratic.Bezier ( 𝕀 Double ) -> 𝕀 Double -> 𝕀 Double +evaluateQuadratic bez t = + -- assert (inf t >= 0 && sup t <= 1) "evaluateCubic: t ⊊ [0,1]" $ -- Requires t ⊆ [0,1] + let inf_bez = fmap inf bez + sup_bez = fmap sup bez + mins = fmap (Quadratic.bezier @( T Double ) inf_bez) + $ inf t :| ( sup t : filter ( inside t ) ( Quadratic.extrema inf_bez ) ) + maxs = fmap (Quadratic.bezier @( T Double ) sup_bez) + $ inf t :| ( sup t : filter ( inside t ) ( Quadratic.extrema sup_bez ) ) + in 𝕀 ( minimum mins ) ( maximum maxs ) diff --git a/src/splines/Math/Interval.hs b/src/splines/Math/Interval.hs index 7920315..a904273 100644 --- a/src/splines/Math/Interval.hs +++ b/src/splines/Math/Interval.hs @@ -12,11 +12,15 @@ module Math.Interval , scaleInterval , 𝕀ℝ , singleton, nonDecreasing + , inside + , aabb ) where -- base import Prelude hiding ( Num(..) ) +import Data.List.NonEmpty + ( NonEmpty(..) ) -- acts import Data.Act @@ -32,20 +36,18 @@ import Data.Group.Generics -- splines import Math.Algebra.Dual -import Math.Algebra.Dual.Internal - ( chainRule1NQ ) import Math.Interval.Internal ( 𝕀(𝕀), inf, sup, scaleInterval ) import Math.Linear ( ℝ(..), T(..) - , RepresentableQ(..) + , RepresentableQ(..), Representable(..) ) import Math.Module import Math.Monomial import Math.Ring -------------------------------------------------------------------------------- --- Interval arithmetic using rounded-hw library. +-- Interval arithmetic. type 𝕀ℝ n = 𝕀 ( ℝ n ) type instance D k ( 𝕀 v ) = D k v @@ -57,6 +59,8 @@ singleton a = 𝕀 a a nonDecreasing :: ( a -> b ) -> 𝕀 a -> 𝕀 b nonDecreasing f ( 𝕀 lo hi ) = 𝕀 ( f lo ) ( f hi ) +inside :: Ord a => 𝕀 a -> a -> Bool +inside ( 𝕀 lo hi ) x = x >= lo && x <= hi deriving via ViaAbelianGroup ( T ( 𝕀 Double ) ) instance Semigroup ( T ( 𝕀 Double ) ) @@ -71,12 +75,19 @@ instance Torsor ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where a --> b = T $ b - a ------------------------------------------------------------------------------- +-- Lattices. + +aabb :: ( Representable r v, Ord r, Functor f ) + => f ( 𝕀 v ) -> ( f ( 𝕀 r ) -> 𝕀 r ) -> 𝕀 v +aabb fv extrema = tabulate \ i -> extrema ( fmap ( `index` i ) fv ) + +-------------------------------------------------------------------------------- instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 0 ) ) where - origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] ) - T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] ) - T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] ) - k *^ T a = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] ) + origin = T $ singleton ℝ0 + _ ^+^ _ = T $ singleton ℝ0 + _ ^-^ _ = T $ singleton ℝ0 + _ *^ _ = T $ singleton ℝ0 instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 1 ) ) where origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] ) diff --git a/src/splines/Math/Ring.hs b/src/splines/Math/Ring.hs index e1d088a..1a8a8bb 100644 --- a/src/splines/Math/Ring.hs +++ b/src/splines/Math/Ring.hs @@ -1,7 +1,8 @@ {-# LANGUAGE ScopedTypeVariables #-} module Math.Ring - ( AbelianGroup(..), Signed(..), Ring(..), Field(..), Transcendental(..) + ( AbelianGroup(..), Signed(..), Ring(..), Field(..) + , Floating(..), Transcendental(..) , ViaPrelude(..), ViaAbelianGroup(..) @@ -11,7 +12,7 @@ module Math.Ring -- base import Prelude ( Num, Fractional ) -import Prelude hiding ( Num(..), Fractional(..) ) +import Prelude hiding ( Num(..), Fractional(..), Floating(..) ) import qualified Prelude import Control.Applicative ( liftA2 ) @@ -60,11 +61,13 @@ class Ring a => Field a where recip x = fromInteger 1 / x x / y = x * recip y +class Field a => Floating a where + sqrt :: a -> a + class Field a => Transcendental a where pi :: a cos :: a -> a sin :: a -> a --- sqrt :: a -> a -------------------------------------------------------------------------------- @@ -126,11 +129,13 @@ instance Fractional a => Field ( ViaPrelude a ) where fromRational = coerce $ Prelude.fromRational @a (/) = coerce $ (Prelude./) @a -instance Floating a => Transcendental ( ViaPrelude a ) where +instance Prelude.Floating a => Floating ( ViaPrelude a ) where + sqrt = coerce $ Prelude.sqrt @a + +instance Prelude.Floating a => Transcendental ( ViaPrelude a ) where pi = coerce $ Prelude.pi @a sin = coerce $ Prelude.sin @a cos = coerce $ Prelude.cos @a --- sqrt = coerce $ Prelude.sqrt @a -------------------------------------------------------------------------------- @@ -160,12 +165,14 @@ deriving via ViaPrelude Float instance AbelianGroup Float deriving via ViaPrelude Float instance Ring Float deriving via ViaPrelude Float instance Signed Float deriving via ViaPrelude Float instance Field Float +deriving via ViaPrelude Float instance Floating Float deriving via ViaPrelude Float instance Transcendental Float deriving via ViaPrelude Double instance AbelianGroup Double deriving via ViaPrelude Double instance Ring Double deriving via ViaPrelude Double instance Signed Double deriving via ViaPrelude Double instance Field Double +deriving via ViaPrelude Double instance Floating Double deriving via ViaPrelude Double instance Transcendental Double -------------------------------------------------------------------------------- diff --git a/src/splines/Math/Roots.hs b/src/splines/Math/Roots.hs index 2f679c7..f713167 100644 --- a/src/splines/Math/Roots.hs +++ b/src/splines/Math/Roots.hs @@ -62,26 +62,26 @@ solveQuadratic -> a -- ^ linear coefficient -> a -- ^ quadratic coefficient -> [ a ] -solveQuadratic a0 a1 a2 - | nearZero a1 && nearZero a2 - = if nearZero a0 +solveQuadratic c b a + | nearZero b && nearZero a + = if nearZero c then [ 0, 0.5, 1 ] -- convention else [] - | nearZero ( a0 * a0 * a2 / ( a1 * a1 ) ) - = [ -a0 / a1 ] + | nearZero ( c * c * a / ( b * b ) ) + = [ -c / b ] | disc < 0 = [] -- non-real solutions | otherwise = let r :: a r = - if a1 >= 0 - then 2 * a0 / ( -a1 - sqrt disc ) - else 0.5 * ( -a1 + sqrt disc ) / a2 - in [ r, -r - a1 / a2 ] + if b >= 0 + then 2 * c / ( -b - sqrt disc ) + else 0.5 * ( -b + sqrt disc ) / a + in [ r, c / ( a * r ) ] where disc :: a - disc = a1 * a1 - 4 * a0 * a2 + disc = b * b - 4 * a * c -------------------------------------------------------------------------------- -- Root finding using Laguerre's method