\documentclass{article}
%include polycode.fmt
\usepackage{fontspec}
\usepackage{amsmath}
\usepackage{unicode-math}
\usepackage{lualatex-math}
\setmainfont{latinmodern-math.otf}
\setmathfont{latinmodern-math.otf}
\usepackage{verbatim}
\author{Sophie Taylor}
\title{haskell-clifford: A Haskell Clifford algebra dynamics library}
\begin{document}
So yeah. This is a Clifford number representation. I will fill out the documentation more fully and stuff once the design has stabilised.
I am basing the design of this on Issac Trotts' geometric algebra library.\cite{hga}
Let us begin. We are going to use the Numeric Prelude because it is (shockingly) nicer for numeric stuff.
\begin{code}
\end{code}
%if False
\begin{code}
\end{code}
%endif
Clifford algebras are a module over a ring. They also support all the usual transcendental functions.
\begin{code}
module Numeric.Clifford.Multivector where
import Numeric.Clifford.Blade
import NumericPrelude hiding (iterate, head, map, tail, reverse, scanl, zipWith, drop, (++), filter, null, length, foldr, foldl1, zip, foldl, concat, (!!), concatMap,any, repeat, replicate, elem, replicate, all)
import Algebra.Absolute
import Algebra.Algebraic
import Algebra.Additive
import Algebra.Ring
import Algebra.OccasionallyScalar
import Algebra.ToInteger
import Algebra.Transcendental
import Algebra.ZeroTestable
import Algebra.Module
import Algebra.Field
import Data.Serialize
import MathObj.Polynomial.Core (progression)
import System.IO
import Data.List.Stream
import Data.Permute (sort, isEven)
import Data.List.Ordered
import Data.Ord
import Data.Maybe
import Numeric.Natural
import qualified Data.Vector as V
import NumericPrelude.Numeric (sum)
import qualified NumericPrelude.Numeric as NPN
import Test.QuickCheck
import Math.Sequence.Converge (convergeBy)
import Control.DeepSeq
import Number.Ratio hiding (scale)
import Algebra.ToRational
import qualified GHC.Num as PNum
import Control.Lens hiding (indices)
import Control.Exception (assert)
import Data.Maybe
import Data.Monoid
import Data.Data
import Data.DeriveTH
import GHC.TypeLits
import Data.Word
import Debug.Trace
\end{code}
A multivector is nothing but a linear combination of basis blades.
\begin{code}
data Multivector (n::Nat) f where
BladeSum :: forall n f . (SingI n, Algebra.Field.C f, Ord f) => { _terms :: [Blade n f]} -> Multivector n f
deriving instance Eq (Multivector n f)
deriving instance Ord (Multivector n f)
deriving instance (Show f) => Show (Multivector n f)
dimension :: forall (n::Nat) f. SingI n => Multivector n f -> Natural
dimension _ = toNatural ((fromIntegral $ fromSing (sing :: Sing n))::Word)
terms :: Lens' (Multivector n f) [Blade n f]
terms = lens _terms (\bladeSum v -> bladeSum {_terms = v})
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
mvTerms m = m^.terms
addLikeTerms' = sumLikeTerms . groupLikeTerms
groupLikeTerms ::Eq f => [Blade n f] -> [[Blade n f]]
groupLikeTerms = groupBy (\a b -> a^.indices == b^.indices)
compensatedSum' :: (Algebra.Additive.C f) => [f] -> f
compensatedSum' xs = kahan zero zero xs where
kahan s _ [] = s
kahan s c (x:xs) =
let y = x c
t = s + y
in kahan t ((ts)y) xs
compensatedRunningSum xs=shanksTransformation . map fst $ scanl kahanSum (zero,zero) xs where
kahanSum (s,c) b = (t,newc) where
y = b c
t = s + y
newc = (t s) y
multiplyOutBlades :: (SingI n, Algebra.Ring.C a) => [Blade n a] -> [Blade n a] -> [Blade n a]
multiplyOutBlades x y = [bladeMul l r | l <-x, r <- y]
multiplyList [] = error "Empty list"
multiplyList l = mvNormalForm $ BladeSum listOfBlades where
expandedBlades a = foldl1 multiplyOutBlades a
listOfBlades = expandedBlades $ map mvTerms l
sumList xs = mvNormalForm $ BladeSum $ concat $ map mvTerms xs
sumLikeTerms :: (Algebra.Field.C f, SingI n) => [[Blade n f]] -> [Blade n f]
sumLikeTerms blades = map (\sameIxs -> map bScale sameIxs & compensatedSum' & (\result -> Blade result ((head sameIxs) & bIndices))) blades
instance (Algebra.Field.C f, SingI n, Ord f) => Data.Monoid.Monoid (Sum (Multivector n f)) where
mempty = Data.Monoid.Sum Algebra.Additive.zero
mappend (Data.Monoid.Sum x) (Data.Monoid.Sum y)= Data.Monoid.Sum (x + y)
mconcat = Data.Monoid.Sum . sumList . map getSum
instance (Algebra.Field.C f, SingI n, Ord f) => Data.Monoid.Monoid (Product (Multivector n f)) where
mempty = Product one
mappend (Product x) (Product y) = Product (x * y)
mconcat = Product . foldl (*) one . map getProduct
e :: (Algebra.Field.C f, Ord f, SingI n) => f -> [Natural] -> Multivector n f
s `e` indices = mvNormalForm $ BladeSum [Blade s indices]
scalar s = s `e` []
instance (Control.DeepSeq.NFData f) => Control.DeepSeq.NFData (Multivector n f)
instance (Control.DeepSeq.NFData f) => Control.DeepSeq.NFData (Blade n f)
instance (Algebra.Field.C f, Ord f, SingI n) => Algebra.Additive.C (Multivector n f) where
a + b = mvNormalForm $ BladeSum (mvTerms a ++ mvTerms b)
a b = mvNormalForm $ BladeSum (mvTerms a ++ map bladeNegate (mvTerms b))
zero = BladeSum [scalarBlade Algebra.Additive.zero]
\end{code}
Now it is time for the Clifford product. :3
\begin{code}
instance (Algebra.Field.C f, Ord f, SingI n) => Algebra.Ring.C (Multivector n f) where
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]
one = scalar Algebra.Ring.one
fromInteger i = scalar $ Algebra.Ring.fromInteger i
a ^ 2 = a * a
a ^ 0 = one
a ^ 1 = a
a ^ n = multiplyList (replicate (NPN.fromInteger n) a)
two = fromInteger 2
mul = (Algebra.Ring.*)
\end{code}
Clifford numbers have a magnitude and absolute value:
\begin{code}
magnitude = sqrt . compensatedSum' . map (\b -> (bScale b)^ 2) . mvTerms
instance (Algebra.Absolute.C f, Algebra.Algebraic.C f, Ord f, SingI n) => Algebra.Absolute.C (Multivector n f) where
abs v = magnitude v `e` []
signum (BladeSum [Blade scale []]) = scalar $ signum scale
signum (BladeSum []) = scalar Algebra.Additive.zero
instance (Algebra.Field.C f, Ord f, SingI n) => Algebra.Module.C f (Multivector n f) where
(*>) s v = v & mvTerms & map (bladeScaleLeft s) & BladeSum
(/) :: (Algebra.Field.C f, Ord f, SingI n) => Multivector n f -> f -> Multivector n f
(/) v d = BladeSum $ map (bladeScaleLeft (NPN.recip d)) $ mvTerms v
(</) n d = Numeric.Clifford.Multivector.inverse d * n
(/>) n d = n * Numeric.Clifford.Multivector.inverse d
(</>) n d = n /> d
integratePoly c x = c : zipWith (Numeric.Clifford.Multivector./) x progression
converge [] = error "converge: empty list"
converge xs = fromMaybe empty (convergeBy checkPeriodic Just xs)
where
empty = error "converge: error in implmentation"
checkPeriodic (a:b:c:_)
| (trace ("Converging at " ++ show a) a) == b = Just a
| a == c = Just a
checkPeriodic _ = Nothing
aitkensAcceleration [] = []
aitkensAcceleration a@(xn:[]) = a
aitkensAcceleration a@(xn:xnp1:[]) = a
aitkensAcceleration a@(xn:xnp1:xnp2:[]) = a
aitkensAcceleration (xn:xnp1:xnp2:xs) | xn == xnp1 = [xnp1]
| xn == xnp2 = [xnp2]
| otherwise = xn ((dxn ^ 2) /> ddxn) : aitkensAcceleration (xnp1:xnp2:xs) where
dxn = sumList [xnp1,negate xn]
ddxn = sumList [xn, (2) * xnp1, xnp2]
shanksTransformation [] = []
shanksTransformation a@(xnm1:[]) = a
shanksTransformation a@(xnm1:xn:[]) = a
shanksTransformation (xnm1:xn:xnp1:xs) | xnm1 == xn = [xn]
| xnm1 == xnp1 = [xnm1]
| denominator == zero = [xnp1]
| otherwise = trace ("Shanks transformation input = " ++ show xn ++ "\nShanks transformation output = " ++ show out) out:shanksTransformation (xn:xnp1:xs) where
out = numerator /> denominator
numerator = sumList [xnp1*xnm1, negate (xn^2)]
denominator = sumList [xnp1, (2)*xn, xnm1]
exp (BladeSum [ Blade s []]) = trace ("scalar exponential of " ++ show s) scalar $ Algebra.Transcendental.exp s
exp x = trace ("Computing exponential of " ++ show x) convergeTerms x where
expMag = Algebra.Transcendental.exp mag
expScaled = converge $ shanksTransformation.shanksTransformation . compensatedRunningSum $ expTerms scaled
convergeTerms terms = converge $ shanksTransformation.shanksTransformation.compensatedRunningSum $ expTerms terms
mag = trace ("In exponential, magnitude is " ++ show ( magnitude x)) magnitude x
scaled = let val = (Numeric.Clifford.Multivector./) x mag in trace ("In exponential, scaled is" ++ show val) val
takeEvery nth xs = case drop (nth1) xs of
(y:ys) -> y : takeEvery nth ys
[] -> []
cosh x = converge $ shanksTransformation . compensatedRunningSum $ takeEvery 2 $ expTerms x
sinh x = converge $ shanksTransformation . compensatedRunningSum $ takeEvery 2 $ tail $ expTerms x
seriesPlusMinus (x:y:rest) = x:Algebra.Additive.negate y: seriesPlusMinus rest
seriesMinusPlus (x:y:rest) = Algebra.Additive.negate x : y : seriesMinusPlus rest
sin x = converge $ shanksTransformation $ compensatedRunningSum $ sinTerms x
sinTerms x = seriesPlusMinus $ takeEvery 2 $ expTerms x
cos x = converge $ shanksTransformation $ compensatedRunningSum (Algebra.Ring.one : cosTerms x)
cosTerms x = seriesMinusPlus $ takeEvery 2 $ tail $ expTerms x
expTerms x = map snd $ iterate (\(n,b) -> (n + 1, (x*b) Numeric.Clifford.Multivector./ fromInteger n )) (1::NPN.Integer,one)
dot a b = mvNormalForm $ BladeSum [x `bDot` y | x <- mvTerms a, y <- mvTerms b]
wedge a b = mvNormalForm $ BladeSum [x `bWedge` y | x <- mvTerms a, y <- mvTerms b]
(∧) = wedge
(⋅) = dot
reverseBlade b = bladeNormalForm $ b & indices %~ reverse
reverseMultivector v = mvNormalForm $ v & terms.traverse%~ reverseBlade
inverse a = assert (a /= zero) $ reverseMultivector a Numeric.Clifford.Multivector./ bScale (head $ mvTerms (a * reverseMultivector a))
recip=Numeric.Clifford.Multivector.inverse
instance (Algebra.Field.C f, Ord f, SingI n) => Algebra.OccasionallyScalar.C f (Multivector n f) where
toScalar = bScale . bladeGetGrade 0 . head . mvTerms
toMaybeScalar (BladeSum [Blade s []]) = Just s
toMaybeScalar (BladeSum []) = Just Algebra.Additive.zero
toMaybeScalar _ = Nothing
fromScalar = scalar
\end{code}
Also, we may as well implement the standard prelude Num interface.
\begin{code}
instance (Algebra.Algebraic.C f, SingI n, Ord f) => PNum.Num (Multivector n f) where
(+) = (Algebra.Additive.+)
() = (Algebra.Additive.-)
(*) = (Algebra.Ring.*)
negate = NPN.negate
abs = scalar . magnitude
fromInteger = Algebra.Ring.fromInteger
signum m = Numeric.Clifford.Multivector.inverse (scalar $ magnitude m) * m
\end{code}
Let's use Newton or Halley iteration to find the principal n-th root :3
\begin{code}
root :: (Show f, Ord f, Algebra.Algebraic.C f, SingI d) => NPN.Integer -> Multivector d f -> Multivector d f
root n (BladeSum [Blade s []]) = scalar $ Algebra.Algebraic.root n s
root n a@(BladeSum _) = converge $ rootIterationsStart n a one
rootIterationsStart ::(Ord f, Show f, Algebra.Algebraic.C f)=> NPN.Integer -> Multivector d f -> Multivector d f -> [Multivector d f]
rootIterationsStart n a@(BladeSum (Blade s [] :xs)) one = rootHalleysIterations n a g where
g = if s >= NPN.zero then one else Algebra.Ring.one `e` [1,2]
rootIterationsStart n a@(BladeSum _) g = rootHalleysIterations n a g
rootNewtonIterations :: (Algebra.Field.C f, Ord f, SingI d) => NPN.Integer -> Multivector d f -> Multivector d f -> [Multivector d f]
rootNewtonIterations n a = iterate xkplus1 where
xkplus1 xk = xk + deltaxk xk
deltaxk xk = oneOverN * ((Numeric.Clifford.Multivector.inverse (xk ^ (n one))* a) xk)
oneOverN = scalar $ NPN.recip $ fromInteger n
rootHalleysIterations :: (Show a, Ord a, Algebra.Algebraic.C a, SingI d) => NPN.Integer -> Multivector d a -> Multivector d a -> [Multivector d a]
rootHalleysIterations n a = halleysMethod f f' f'' where
f x = a (x^n)
f' x = fromInteger (n) * (x^(n1))
f'' x = fromInteger ((n*(n1))) * (x^(n2))
pow a p = (a ^ up) Numeric.Clifford.Multivector./> Numeric.Clifford.Multivector.root down a where
ratio = toRational p
up = numerator ratio
down = denominator ratio
halleysMethod :: (Show a, Ord a, Algebra.Algebraic.C a, SingI d) => (Multivector d a -> Multivector d a) -> (Multivector d a -> Multivector d a) -> (Multivector d a -> Multivector d a) -> Multivector d a -> [Multivector d 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)
fx = f x
dfx = f' x
ddfx = f'' x
secantMethod f x0 x1 = update x1 x0 where
update xm1 xm2 | xm1 == xm2 = [xm1]
| otherwise = if x == xm1 then [x] else x : update x xm1 where
x = xm1 f xm1 * (xm1xm2) * Numeric.Clifford.Multivector.inverse (f xm1 f xm2)
\end{code}
Now let's try logarithms by fixed point iteration. It's gonna be slow, but whatever!
\begin{code}
normalised a = a * (scalar $ NPN.recip $ magnitude a)
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 = converge $ halleysMethod f f' f'' (one `e` [1,2]) where
f x = a Numeric.Clifford.Multivector.exp x
f' x = NPN.negate $ Numeric.Clifford.Multivector.exp x
f'' = f'
\end{code}
Now let's do (slow as fuck probably) numerical integration! :D~! Since this is gonna be used for physical applications, it's we're gonna start off with a Hamiltonian structure and then a symplectic integrator.
\begin{code}
\end{code}
\bibliographystyle{IEEEtran}
\bibliography{biblio.bib}
\end{document}