clifford 0.1.0.8 → 0.1.0.9
raw patch · 9 files changed
+130/−45 lines, 9 files
Files
- bench/benchmarks.hs +1/−2
- changelog.md +1/−0
- clifford.cabal +1/−1
- src/Numeric/Clifford/Blade.lhs +20/−11
- src/Numeric/Clifford/Internal.hs +10/−0
- src/Numeric/Clifford/LinearOperators.lhs +1/−1
- src/Numeric/Clifford/Multivector.lhs +82/−29
- src/Numeric/Clifford/NumericIntegration.lhs +13/−0
- test/Numeric/Clifford/MultivectorSpec.lhs +1/−1
bench/benchmarks.hs view
@@ -14,8 +14,7 @@ import NumericPrelude hiding (iterate, last, map, take, log) import Prelude hiding (iterate, last, map, negate, take,log, (*), (+))- -type STVector = Multivector 3 1 Double+ scalar2 = scalar (2.0::NumericPrelude.Double) :: STVector ij2 = (2.0::NumericPrelude.Double) `e` [1,2] :: STVector ik3 = (3::NumericPrelude.Double) `e` [1,3] :: STVector
changelog.md view
@@ -1,4 +1,5 @@ -*-change-log-*-+ 0.1.0.9 Inlined/specialised a bunch of function, hueg speed increase 0.1.0.8 Implemented algebraic/transcendental typeclasses 0.1.0.7 Adding basic linear operators; made multivector a field 0.1.0.6 Memoising the blade index comparision function for a 20% speed increase
clifford.cabal view
@@ -10,7 +10,7 @@ -- PVP summary: +-+------- breaking API changes -- | | +----- non-breaking API additions -- | | | +--- code changes with no API change-version: 0.1.0.8+version: 0.1.0.9 -- A short (one-line) description of the package. synopsis: A Clifford algebra library
src/Numeric/Clifford/Blade.lhs view
@@ -66,13 +66,18 @@ data Blade (p :: Nat) (q :: Nat) f where Blade :: forall p q f . (SingI p, SingI q, Algebra.Field.C f) => {_scale :: f, _indices :: [Natural]} -> Blade p q f +type STBlade = Blade 3 1 Double+type E3Blade = Blade 3 0 Double scale :: Lens' (Blade p q f) f scale = lens _scale (\blade v -> blade {_scale = v}) indices :: Lens' (Blade p q f) [Natural] indices = lens _indices (\blade v -> blade {_indices = v}) dimension :: forall (p::Nat) (q::Nat) f. (SingI p, SingI q) => Blade p q f -> (Natural,Natural) dimension _ = (toNatural ((GHC.Real.fromIntegral $ fromSing (sing :: Sing p))::Word),toNatural((GHC.Real.fromIntegral $ fromSing (sing :: Sing q))::Word))++bScale :: Blade p q f -> f bScale b = b^.scale+bIndices :: Blade p q f -> [Natural] bIndices b = b^.indices instance (Control.DeepSeq.NFData f) => Control.DeepSeq.NFData (Blade p q f) instance(Show f) => Show (Blade p q f) where@@ -99,11 +104,15 @@ zeroBlade :: (Algebra.Field.C f, SingI p, SingI q) => Blade p q f zeroBlade = scalarBlade Algebra.Additive.zero +bladeNonZero :: (Algebra.Additive.C f, Eq f) => Blade p q f -> Bool bladeNonZero b = b^.scale /= Algebra.Additive.zero +bladeNegate :: (Algebra.Additive.C f) => Blade p q f -> Blade p q f bladeNegate b = b&scale%~negate --Blade (Algebra.Additive.negate$ b^.scale) (b^.indices) +bladeScaleLeft :: f -> Blade p q f -> Blade p q f bladeScaleLeft s (Blade f ind) = Blade (s * f) ind+bladeScaleRight :: f -> Blade p q f -> Blade p q f bladeScaleRight s (Blade f ind) = Blade (f * s) ind \end{code} @@ -118,6 +127,9 @@ \begin{code} +{-#INLINE bladeNormalForm#-}+{-#SPECIALISE INLINE bladeNormalForm::E3Blade -> E3Blade #-}+{-#SPECIALISE INLINE bladeNormalForm :: STBlade -> STBlade #-} bladeNormalForm :: forall (p::Nat) (q::Nat) f. Blade p q f -> Blade p q f bladeNormalForm (Blade scale indices) = result where@@ -130,6 +142,7 @@ scale' = if doNotNegate then scale else negate scale (newIndices, doNotNegate) = sortIndices (indices,q') +sortIndices :: ([Natural],Integer) -> ([Natural],Bool) sortIndices = memo sortIndices' where sortIndices' :: ([Natural],Integer) -> ([Natural],Bool) sortIndices' (indices,q') = (uniqueSorted, doNotNegate) where@@ -167,6 +180,9 @@ First up for operations: Blade multiplication. This is no more than assembling orthogonal vectors into k-vectors. \begin{code}+{-#INLINE bladeMul #-}+{-#SPECIALISE INLINE bladeMul :: STBlade -> STBlade -> STBlade #-}+{-#SPECIALISE INLINE bladeMul :: E3Blade -> E3Blade -> E3Blade #-} bladeMul :: Blade p q f -> Blade p q f-> Blade p q f bladeMul x@(Blade _ _) y@(Blade _ _)= bladeNormalForm $ Blade (bScale x Algebra.Ring.* bScale y) (bIndices x ++ bIndices y) multiplyBladeList :: (SingI p, SingI q, Algebra.Field.C f) => [Blade p q f] -> Blade p q f@@ -200,6 +216,7 @@ k = Algebra.Absolute.abs $ grade x - grade y xy = bladeMul x y +propBladeDotAssociative :: (Algebra.Additive.C f, Eq f) => Blade p q f -> Blade p q f -> Blade p q f -> Bool propBladeDotAssociative = Algebra.Laws.associative bDot \end{code}@@ -212,17 +229,14 @@ instance (Algebra.Additive.C f, Ord f) => Ord (Blade p q f) where compare a b | bIndices a == bIndices b = compare (bScale a) (bScale b) | otherwise = compareIndices (bIndices a) (bIndices b)++compareIndices :: [Natural] -> [Natural] -> Ordering compareIndices = memo compareIndices' where compareIndices' a b = case compare (length a) (length b) of LT -> LT GT -> GT EQ -> compare a b -instance Arbitrary Natural where- arbitrary = sized $ \n ->- let n' = NPN.abs n in- fmap (toNatural . (\x -> (GHC.Real.fromIntegral x)::Word)) (choose (0, n'))- shrink = shrinkIntegral instance (SingI p, SingI q, Algebra.Field.C f, Arbitrary f) => Arbitrary (Blade p q f) where arbitrary = do@@ -247,13 +261,8 @@ \begin{code} -{- Note: Figure out what this is meant to be lol-skewcommutative op x y = x `op` y == (bladeScaleLeft (fromInteger (-1))$ y `op` x) -propAnticommutativeMultiplication :: (Eq f,Algebra.Ring.C f, Algebra.Additive.C f) => Blade f -> Blade f -> Bool-propAnticommutativeMultiplication = anticommutative bladeMul--}-propCommutativeAddition = commutative (+)+--propCommutativeAddition = commutative (+) \end{code} \bibliographystyle{IEEEtran} \bibliography{biblio.bib}
src/Numeric/Clifford/Internal.hs view
@@ -7,6 +7,8 @@ import Data.List.Stream import Control.Arrow import Data.Bits+import Test.QuickCheck+import Data.Word import qualified Debug.Trace as DebugTrace #ifdef DEBUG myTrace = DebugTrace.trace@@ -18,6 +20,14 @@ trie f = NaturalTrie (trie (f . unbitsZ)) untrie (NaturalTrie t) = untrie t . bitsZ enumerate (NaturalTrie t) = enum' unbitsZ t+++instance Arbitrary Natural where+ arbitrary = sized $ \n ->+ let n' = abs n in+ fmap (toNatural . (\x -> (fromIntegral x)::Word)) (choose (0, n'))+ shrink = shrinkIntegral+ unbitsZ :: (Prelude.Num n, Bits n) => (Bool,[Bool]) -> n
src/Numeric/Clifford/LinearOperators.lhs view
@@ -19,7 +19,7 @@ makeReflectionOperator u = reflect u rotate spinor x = (reverseMultivector spinor) * x * spinor-rotatePlaneAngle plane angle = rotate (exp ((normalised plane) * (angle/2)))+rotatePlaneAngle plane angle = rotate (exp (((fst.normalised) plane) * (angle/2))) makeRotationOperator :: LinearOperatorCreator p q f makeRotationOperator u = rotate u
src/Numeric/Clifford/Multivector.lhs view
@@ -88,6 +88,9 @@ data Multivector (p::Nat) (q::Nat) f where BladeSum :: forall p q f . (Ord f, Algebra.Field.C f, SingI p, SingI q) => { _terms :: [Blade p q f]} -> Multivector p q f +type STVector = Multivector 3 1 Double+type E3Vector = Multivector 3 0 Double+ instance (SingI p, SingI q, Algebra.Field.C f, Arbitrary f, Ord f) => Arbitrary (Multivector p q f) where arbitrary = mvNormalForm <$> BladeSum <$> (vector d) where p' = (fromSing (sing :: Sing p)) :: Integer@@ -107,18 +110,23 @@ terms :: Lens' (Multivector p q f) [Blade p q f] terms = lens _terms (\bladeSum v -> bladeSum {_terms = v}) +{-# INLINE mvNormalForm #-} mvNormalForm (BladeSum terms) = BladeSum $ if null resultant then [scalarBlade Algebra.Additive.zero] else resultant where resultant = filter bladeNonZero $ addLikeTerms' $ Data.List.Ordered.sortBy compare $ map bladeNormalForm $ terms+{-#INLINE mvTerms #-} mvTerms m = m^.terms +{-# INLINE addLikeTerms' #-} addLikeTerms' = sumLikeTerms . groupLikeTerms +{-# INLINE groupLikeTerms #-} groupLikeTerms ::Eq f => [Blade p q f] -> [[Blade p q f]] groupLikeTerms = groupBy (\a b -> a^.indices == b^.indices) compareTol :: (Algebra.Algebraic.C f, Algebra.Absolute.C f, Ord f, SingI p, SingI q) => Multivector p q f -> Multivector p q f -> f -> Bool compareTol x y tol = ((NPN.abs $ magnitude (x-y) ) <= tol) +{-#INLINE compensatedSum' #-} compensatedSum' :: (Algebra.Additive.C f) => [f] -> f compensatedSum' xs = kahan zero zero xs where kahan s _ [] = s@@ -128,7 +136,10 @@ in kahan t ((t-s)-y) xs --use this to sum taylor series et al with converge---compensatedRunningSum :: (Algebra.Additive.C f) => [f] -> [f]+{-#INLINE compensatedRunningSum#-}+{-#SPECIALISE INLINE compensatedRunningSum :: [STVector] -> [STVector] #-}+{-#SPECIALISE INLINE compensatedRunningSum :: [E3Vector] -> [E3Vector] #-}+compensatedRunningSum :: (Algebra.Algebraic.C f, Ord f, SingI p, SingI q, Show f) => [Multivector p q f] -> [Multivector p q f] compensatedRunningSum xs=shanksTransformation . map fst $ scanl kahanSum (zero,zero) xs where kahanSum (s,c) b = (t,newc) where y = b - c@@ -164,8 +175,12 @@ --things to test: is 1. adding blades into a map based on indices 2. adding errything together 3. sort results quicker than -- 1. sorting by indices 2. groupBy-ing on indices 3. adding the lists of identical indices +{-#INLINE sumList #-} sumList xs = mvNormalForm $ BladeSum $ concat $ map mvTerms xs +{-#INLINE sumLikeTerms #-}+{-#SPECIALISE INLINE sumLikeTerms :: [[STBlade]] -> [STBlade] #-}+{-#SPECIALISE INLINE sumLikeTerms :: [[E3Blade]] -> [E3Blade] #-} sumLikeTerms :: (Algebra.Field.C f, SingI p, SingI q) => [[Blade p q f]] -> [Blade p q f] sumLikeTerms blades = map (\sameIxs -> map bScale sameIxs & compensatedSum' & (\result -> Blade result ((head sameIxs) & bIndices))) blades @@ -181,15 +196,29 @@ mconcat = Product . foldl (*) one . map getProduct --Constructs a multivector from a scaled blade.+{-#INLINE e#-} e :: (Algebra.Field.C f, Ord f, SingI p, SingI q) => f -> [Natural] -> Multivector p q f s `e` indices = mvNormalForm $ BladeSum [Blade s indices]+{-#INLINE scalar#-} scalar s = s `e` [] instance (Control.DeepSeq.NFData f) => Control.DeepSeq.NFData (Multivector p q f) ++{-{-# RULES+ "turn multiple additions into sumList" forall (f::Algebra.Field.C) (a::Multivector p q f) b c . (+) a ((+) b c) = sumList [a,b,c]+ #-}-}+{-#RULES+ "sumList[..] + a = sumList [..,a]" forall (a::Multivector (p::Nat) (q::Nat) (Algebra.Field.C f)) xs. (+) (sumList xs) a = sumList (a:xs)+ #-}+{-# RULES+ "a+ sumList[..] = sumList [..,a]" forall (a::Multivector p q (Algebra.Field.C f)) xs. (+) a (sumList xs) = sumList (a:xs)+ #-} instance (Algebra.Field.C f, Ord f, SingI p, SingI q) => Algebra.Additive.C (Multivector p q f) where+ {-#INLINE (+)#-} a + b = mvNormalForm $ BladeSum (mvTerms a ++ mvTerms b)+ {-#INLINE (-)#-} a - b = mvNormalForm $ BladeSum (mvTerms a ++ map bladeNegate (mvTerms b)) zero = BladeSum [scalarBlade Algebra.Additive.zero] @@ -201,6 +230,7 @@ \begin{code} instance (Algebra.Field.C f, Ord f, SingI p, SingI q) => Algebra.Ring.C (Multivector p q f) where+ {-#INLINE (*)#-} BladeSum [Blade s []] * b = BladeSum $ map (bladeScaleLeft s) $ mvTerms b a * BladeSum [Blade s []] = BladeSum $ map (bladeScaleRight s) $ mvTerms a a * b = mvNormalForm $ BladeSum [bladeMul x y | x <- mvTerms a, y <- mvTerms b]@@ -213,6 +243,7 @@ --a ^ n --n < 0 = Clifford.recip $ a ^ (negate n) a ^ n = multiplyList (replicate (NPN.fromInteger n) a) + two = fromInteger 2 mul = (Algebra.Ring.*) @@ -228,7 +259,11 @@ \begin{code} ---magnitude :: (Algebra.Algebraic.C f) => Multivector f -> f++{-# INLINE magnitude #-}+{-# SPECIALISE INLINE magnitude:: Multivector 3 1 Double -> Double #-}+{-# SPECIALISE INLINE magnitude:: Multivector 3 0 Double -> Double #-}+magnitude :: (Algebra.Algebraic.C f) => Multivector p q f -> f magnitude = sqrt . compensatedSum' . map (\b -> (bScale b)^ 2) . mvTerms instance (Algebra.Absolute.C f, Algebra.Algebraic.C f, Ord f, SingI p, SingI q) => Algebra.Absolute.C (Multivector p q f) where@@ -244,8 +279,9 @@ --(/) :: (Algebra.Field.C f, Ord f, SingI p, SingI q) => Multivector p q f -> f -> Multivector p q f --(/) v d = BladeSum $ map (bladeScaleLeft (NPN.recip d)) $ mvTerms v --Algebra.Field.recip d *> v-+{-#INLINE (</)#-} (</) n d = Numeric.Clifford.Multivector.inverse d * n+{-#INLINE (/>)#-} (/>) n d = n * Numeric.Clifford.Multivector.inverse d (</>) n d = n /> d @@ -257,7 +293,7 @@ divideRight v s = scaleRight v (recip s) --integratePoly c x = c : zipWith (Numeric.Clifford.Multivector./) x progression ---converge :: (Eq f, Show f) => [f] -> f+{-# INLINE converge#-} converge [] = error "converge: empty list" converge xs = fromMaybe empty (convergeBy checkPeriodic Just xs) where@@ -279,6 +315,10 @@ dxn = sumList [xnp1,negate xn] ddxn = sumList [xn, (-2) * xnp1, xnp2] +{-# INLINABLE shanksTransformation #-}+{-#SPECIALISE shanksTransformation :: [Multivector 3 0 Double] -> [Multivector 3 0 Double] #-}+{-#SPECIALISE shanksTransformation :: [Multivector 3 1 Double] -> [Multivector 3 1 Double] #-}+shanksTransformation :: (Algebra.Algebraic.C f, Ord f, Show f, SingI p, SingI q) => [Multivector p q f] -> [Multivector p q f] shanksTransformation [] = [] shanksTransformation a@(xnm1:[]) = a shanksTransformation a@(xnm1:xn:[]) = a@@ -291,13 +331,9 @@ denominator = sumList [xnp1, (-2)*xn, xnm1] ---exp ::(Ord f, Show f, Algebra.Transcendental.C f)=> Multivector f -> Multivector f ----+{-# INLINABLE takeEvery #-} takeEvery nth xs = case drop (nth-1) xs of (y:ys) -> y : takeEvery nth ys [] -> []@@ -311,34 +347,45 @@ -+{-#INLINE expTerms#-}+{-# SPECIALISE INLINE expTerms :: STVector -> [STVector]#-}+{-# SPECIALISE INLINE expTerms :: E3Vector -> [E3Vector]#-} expTerms :: (Algebra.Algebraic.C f, SingI p, SingI q, Ord f) => Multivector p q f -> [Multivector p q f] expTerms x = map snd $ iterate (\(n,b) -> (n + 1, (recip $ fromInteger n ) `scaleLeft` (x*b) )) (1::NPN.Integer,one) instance (Algebra.Transcendental.C f, Ord f, SingI p, SingI q, Show f) => Algebra.Transcendental.C (Multivector p q f) where pi = scalar pi- exp (BladeSum [ Blade s []]) = myTrace ("scalar exponential of " ++ show s) scalar $ Algebra.Transcendental.exp s+ {-#INLINABLE exp#-}+ {-# SPECIALISE INLINE exp :: STVector -> STVector #-}+ {-# SPECIALISE INLINE exp :: E3Vector -> E3Vector #-}+ exp (BladeSum [ Blade s []]) = myTrace ("scalar exponential of " ++ show s) scalar $ exp s exp x = myTrace ("Computing exponential of " ++ show x) convergeTerms x where --(expMag ^ expScaled) where- expMag = Algebra.Transcendental.exp mag+ expMag = exp mag expScaled = converge $ shanksTransformation.shanksTransformation . compensatedRunningSum $ expTerms scaled convergeTerms terms = converge $ shanksTransformation.shanksTransformation.compensatedRunningSum $ expTerms terms mag = myTrace ("In exponential, magnitude is " ++ show ( magnitude x)) magnitude x scaled = let val = (recip mag) *> x in myTrace ("In exponential, scaled is" ++ show val) val-+ {-#INLINE log#-}+ {-# SPECIALISE INLINE log :: STVector -> STVector #-}+ {-# SPECIALISE INLINE log :: E3Vector -> E3Vector #-} log (BladeSum [Blade s []]) = scalar $ NPN.log s- log a = scalar (NPN.log mag) + log' scaled where- scaled = normalised a- mag = magnitude a+ log a = scalar (log mag) + log' scaled where+ (scaled,mag) = normalised a log' a = converge $ halleysMethod f f' f'' (one `e` [1,2]) where+ {-#INLINABLE f#-} f x = a - exp x+ {-#INLINABLE f'#-} f' x = NPN.negate $ exp x+ {-#INLINABLE f''#-} f'' = f'-+ sin (BladeSum [Blade s []]) = scalar $ sin s sin x = converge $ shanksTransformation $ compensatedRunningSum $ sinTerms x where sinTerms x = seriesPlusMinus $ takeEvery 2 $ expTerms x+ cos (BladeSum [Blade s []]) = scalar $ cos s cos x = converge $ shanksTransformation $ compensatedRunningSum (one : cosTerms x) where cosTerms x = seriesMinusPlus $ takeEvery 2 $ tail $ expTerms x-+ + atan (BladeSum [Blade s []]) = scalar $ atan s atan z = (z/onePlusZSquared) * (one + (converge $ shanksTransformation $ compensatedRunningSum $ map lambda [1..])) where lambda :: Integer -> Multivector p q f lambda n = multiplyList1 $ map innerFraction [1..n]@@ -357,10 +404,12 @@ (∧) = wedge :: Multivector p q f -> Multivector p q f -> Multivector p q f (⋅) = dot :: Multivector p q f -> Multivector p q f -> Multivector p q f +{-# INLINE reverseBlade #-} reverseBlade b = bladeNormalForm $ b & indices %~ reverse +{-# INLINE reverseMultivector #-} reverseMultivector v = mvNormalForm $ v & terms.traverse%~ reverseBlade -+{-#INLINE inverse#-} inverse a@(BladeSum _) = assert (a /= zero) $ (recip scalarComponent) *> (reverseMultivector a) where scalarComponent = bScale (head $ mvTerms (a * reverseMultivector a)) @@ -398,14 +447,14 @@ root 0 _ = error "Cannot take 0th root" root _ (BladeSum []) = error "Empty bladesum" root _ (BladeSum [Blade zero []]) = error "Cannot compute a root of zero"- root n (BladeSum [Blade s []]) = scalar $ Algebra.Algebraic.root n s+ root n (BladeSum [Blade s []]) = scalar $ root n s root n a@(BladeSum _) = converge $ rootIterationsStart n a g where g = if q' <= 1 then one`e`[q',succ q'] else one + one `e` [0,1] (p',q') = signature a rootIterationsStart ::(Ord f, Show f, Algebra.Algebraic.C f)=> NPN.Integer -> Multivector p q f -> Multivector p q f -> [Multivector p q f]-rootIterationsStart n a@(BladeSum (Blade s [] :xs)) one = rootHalleysIterations n a g where- g = if s >= NPN.zero || q' == 1 then one else Algebra.Ring.one `e` [0,1] +rootIterationsStart n a@(BladeSum (Blade s [] :_)) one = rootHalleysIterations n a g where+ g = if s >= NPN.zero || q' == 1 then one else one `e` [0,1] (p',q') = signature a rootIterationsStart n a@(BladeSum _) g = rootHalleysIterations n a g@@ -414,7 +463,7 @@ rootNewtonIterations :: (Algebra.Field.C f, Ord f, SingI p, SingI q) => NPN.Integer -> Multivector p q f -> Multivector p q f -> [Multivector p q f] rootNewtonIterations n a = iterate xkplus1 where xkplus1 xk = xk + deltaxk xk- deltaxk xk = oneOverN * ((Numeric.Clifford.Multivector.inverse (xk ^ (n - one))* a) - xk)+ deltaxk xk = oneOverN * ((inverse (xk ^ (n - one))* a) - xk) oneOverN = scalar $ NPN.recip $ fromInteger n rootHalleysIterations :: (Show a, Ord a, Algebra.Algebraic.C a, SingI p, SingI q) => NPN.Integer -> Multivector p q a -> Multivector p q a -> [Multivector p q a]@@ -429,11 +478,12 @@ up = numerator ratio down = denominator ratio-} +{-#INLINE halleysMethod #-} halleysMethod :: (Show a, Ord a, Algebra.Algebraic.C a, SingI p, SingI q) => (Multivector p q a -> Multivector p q a) -> (Multivector p q a -> Multivector p q a) -> (Multivector p q a -> Multivector p q a) -> Multivector p q a -> [Multivector p q a] halleysMethod f f' f'' = iterate update where- update x = x - (numerator x * Numeric.Clifford.Multivector.inverse (denominator x)) where- numerator x = multiplyList [2, fx, dfx]- denominator x = multiplyList [2, dfx, dfx] - (fx * ddfx)+ update x = x - (numerator x * inverse (denominator x) ) where+ numerator x= multiplyList [2, fx, dfx]+ denominator x= multiplyList [2, dfx, dfx] - (fx * ddfx) fx = f x dfx = f' x ddfx = f'' x@@ -450,9 +500,12 @@ Now let's try logarithms by fixed point iteration. It's gonna be slow, but whatever! \begin{code}--normalised :: (Ord f, Algebra.Algebraic.C f, SingI p, SingI q) => Multivector p q f -> Multivector p q f-normalised a = a `scaleRight` ( NPN.recip $ magnitude a)+{-#INLINE normalised#-}+{-#SPECIALISE INLINE normalised :: STVector -> (STVector, Double) #-}+{-#SPECIALISE INLINE normalised :: E3Vector -> (E3Vector, Double) #-}+normalised :: (Ord f, Algebra.Algebraic.C f, SingI p, SingI q) => Multivector p q f -> (Multivector p q f,f)+normalised a = (a `scaleRight` ( recip $ mag),mag) where+ mag = magnitude a \end{code}
src/Numeric/Clifford/NumericIntegration.lhs view
@@ -134,6 +134,10 @@ $( derive makeIs ''RKAttribute) +{-#SPECIALISE genericRKMethod :: ButcherTableau Double -> [RKAttribute Double stateType] -> RKStepper 3 0 Double stateType#-}+{-#SPECIALISE genericRKMethod :: ButcherTableau Double -> [RKAttribute Double [E3Vector]] -> RKStepper 3 0 Double [E3Vector]#-}+{-#SPECIALISE genericRKMethod :: ButcherTableau Double -> [RKAttribute Double stateType] -> RKStepper 3 1 Double stateType#-}+{-#SPECIALISE genericRKMethod :: ButcherTableau Double -> [RKAttribute Double [STVector]] -> RKStepper 3 1 Double [STVector]#-} genericRKMethod :: forall (p::Nat) (q::Nat) t stateType . ( Ord t, Show t, Algebra.Module.C t (Multivector p q t),Algebra.Absolute.C t, Algebra.Algebraic.C t, SingI p, SingI q) => ButcherTableau t -> [RKAttribute t stateType] -> RKStepper p q t stateType@@ -144,6 +148,7 @@ c n = l !! (n-1) where l = _tableauC tableau a :: Int -> [t]+ {-#INLINE a#-} a n = (l !! (n-1)) & filter (/= zero) where l = _tableauA tableau b :: Int -> t@@ -164,6 +169,11 @@ Just (AdaptiveStepSize sigma) -> sigma Nothing -> (\_ _ -> one) + {-#INLINE rkMethodImplicitFixedPoint#-}+-- {-#SPECIALISE rkMethodImplicitFixedPoint :: RKStepper 3 0 Double stateType #-}+-- {-#SPECIALISE rkMethodImplicitFixedPoint :: RKStepper 3 0 Double [E3Vector] #-}+-- {-#SPECIALISE rkMethodImplicitFixedPoint :: RKStepper 3 1 Double stateType #-}+-- {-#SPECIALISE rkMethodImplicitFixedPoint :: RKStepper 3 1 Double [STVector] #-} rkMethodImplicitFixedPoint :: RKStepper p q t stateType rkMethodImplicitFixedPoint h f project unproject (time, state) = (time + (stepSizeAdapter time state)*h*(c s), newState) where@@ -180,8 +190,10 @@ guessTime = time + h' zkp1 :: NPN.Int -> [Multivector p q t] -> [Multivector p q t] zkp1 i zk = map (h*>) (sumOfJs i zk) where+ {-#INLINE sumOfJs#-} sumOfJs :: Int -> [Multivector p q t] -> [Multivector p q t] sumOfJs i zk = sumListOfLists $ map (scaledByAij zk) (a i) where + {-# INLINE scaledByAij #-} scaledByAij :: [Multivector p q t] -> t -> [Multivector p q t] scaledByAij guess a = map (a*>) $ evalDerivatives guessTime $ elementAdd state' guess @@ -189,6 +201,7 @@ newState :: stateType newState = project $ elementAdd state' (assert (not $ null dy) dy) dy = sumListOfLists [map ((b i) *>) (zi i) | i <- [1..s]] :: [Multivector p q t]+ {-#INLINE evalDerivatives #-} evalDerivatives :: t -> [Multivector p q t] -> [Multivector p q t] evalDerivatives time stateAtTime= unproject $ (f time) $ project stateAtTime
test/Numeric/Clifford/MultivectorSpec.lhs view
@@ -13,7 +13,7 @@ main :: IO () main = hspec spec -type STVector = Multivector 3 1 Double+ spec :: Spec spec = do let i = 1.0 `e` [1,2] :: STVector