diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,20 @@ +Copyright (c) 2016-2017 Matti A. Eskelinen + +Permission is hereby granted, free of charge, to any person obtaining +a copy of this software and associated documentation files (the +"Software"), to deal in the Software without restriction, including +without limitation the rights to use, copy, modify, merge, publish, +distribute, sublicense, and/or sell copies of the Software, and to +permit persons to whom the Software is furnished to do so, subject to +the following conditions: + +The above copyright notice and this permission notice shall be included +in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY +CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, +TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE +SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/benchmarks/bench-clifProduct.hs b/benchmarks/bench-clifProduct.hs new file mode 100644 --- /dev/null +++ b/benchmarks/bench-clifProduct.hs @@ -0,0 +1,30 @@ +{- + - A simple benchmark for the geometric product. +-} +{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, BangPatterns, RankNTypes #-} +module Main where +import Clif +import Data.List + +import Data.Time.Clock + +alphabet = [1 .. 20] :: [Int] +blades = subsequences alphabet + +eBlades = map ((1 *:) . map E) blades :: [Clif (Euclidean Int) Int] +lBlades = map ((1 *:) . map S) blades :: [Clif (Lorentzian Int) Int] + +prods bs = zipWith (*) bs bs + +main :: IO () +main = do + putStrLn $ "Calculating the sum of " ++ show (length blades) ++ " geometric products took ..." + start1 <- getCurrentTime + let !s1 = sum (prods eBlades) + end1 <- getCurrentTime + putStrLn $ show (diffUTCTime end1 start1) ++ " with Euclidean Int basis" + + start2 <- getCurrentTime + let !s2 = sum (prods lBlades) + end2 <- getCurrentTime + putStrLn $ show (diffUTCTime end2 start2)++ " with Lorentzian Int basis" diff --git a/clif.cabal b/clif.cabal new file mode 100644 --- /dev/null +++ b/clif.cabal @@ -0,0 +1,65 @@ +name: clif +version: 0.1.0.0 +synopsis: A Clifford algebra number type for Haskell +description: Clif is a library for symbolic and numeric computations on Clifford algebras using Haskell. It is general enough to handle finite and infinite-dimensional Clifford algebras arising from arbitrary bilinear forms, within limitations of computability. To get started, read "Clif.Tutorial". +license: MIT +license-file: LICENSE +author: Matti Eskelinen +maintainer: matti.a.eskelinen@gmail.com +copyright: (c) 2016-2017 Matti Eskelinen +category: Math, Algebra +build-type: Simple +-- extra-source-files: +cabal-version: >=1.10 + +source-repository head + type: git + location: git@github.com:maaleske/clif.git + +source-repository this + type: git + location: git@github.com:maaleske/clif.git + tag: 0.1.0.0 + +library + exposed-modules: Clif + , Clif.Algebra + , Clif.Basis + , Clif.Internal + , Clif.Arbitrary + , Clif.Tutorial + other-modules: + ghc-options: -Wall + other-extensions: FlexibleInstances + , MultiParamTypeClasses + , GeneralizedNewtypeDeriving + , DeriveGeneric + build-depends: base >=4.8 && <4.10 + , QuickCheck >=2.8 && <2.9 + , containers >=0.5 && <0.6 + hs-source-dirs: src + default-language: Haskell2010 + +test-suite Tests + hs-source-dirs: tests + main-is: tests.hs + ghc-options: -Wall + type: exitcode-stdio-1.0 + build-depends: clif + , base >=4.8 && <4.10 + , containers >=0.5 && <0.6 + , QuickCheck >=2.8 && <2.9 + , tasty + , tasty-th + , tasty-quickcheck + default-language: Haskell2010 + +Benchmark bench-clifProduct + hs-source-dirs: benchmarks + main-is: bench-clifProduct.hs + ghc-options: -Wall -O2 + type: exitcode-stdio-1.0 + build-depends: clif + , base>=4.8 && <4.10 + , time>=1.5 && <1.8 + default-language: Haskell2010 diff --git a/src/Clif.hs b/src/Clif.hs new file mode 100644 --- /dev/null +++ b/src/Clif.hs @@ -0,0 +1,21 @@ +{-| +Module : Clif +Copyright : (c) Matti A. Eskelinen 2016-2017 +License : MIT +Maintainer : matti.a.eskelinen@gmail.com +Stability : experimental +Portability : POSIX + +This helper module exports the main modules needed to construct and operate on Clifford algebra values. If you just want to use the library for computations, this should be all you need. To get started, read "Clif.Tutorial". + +-} +module Clif + ( + -- * Construction of a basis + module Clif.Basis + -- * Algebraic operations + , module Clif.Algebra + + ) where +import Clif.Basis +import Clif.Algebra diff --git a/src/Clif/Algebra.hs b/src/Clif/Algebra.hs new file mode 100644 --- /dev/null +++ b/src/Clif/Algebra.hs @@ -0,0 +1,99 @@ +{-| +Module : Clif.Algebra +Copyright : (c) Matti A. Eskelinen, 2016-2017 +License : MIT +Maintainer : matti.a.eskelinen@gmail.com +Stability : experimental +Portability : POSIX + +This module provides the type 'Clif' for representing the elements of a Clifford algebra along with some standard operations. + +See /The inner products of geometric algebra/ by Leo Dorst for a concise explanation of the different inner products. +-} +{-# LANGUAGE + FlexibleInstances, + MultiParamTypeClasses + #-} +module Clif.Algebra + ( + -- * The @Clif@ type + Clif + + -- * Constructing and deconstructing @Clifs@ + , blade, (*:) + , vec + , fromList + + -- * Geometric algebra operations + , grade, rev + + -- * Outer product + , wedge, (/\) + + -- * Inner products + , (<\), (/>), (.|.), (<.>) + , lContract, rContract, scalarProd, dot, hestenes + + -- * Hodge duality + , hodge + + -- * Projections + , proj + + ) where +import Clif.Basis +import Clif.Internal + +infixl 8 <\, />, .|., <.>, /\ +(<\), (/>), (.|.), (<.>), (/\) :: (Eq a, Basis b a) => Clif b a -> Clif b a -> Clif b a +-- | Infix synonym for 'lContract' +(<\) = lContract +-- | Infix synonym for 'rContract' +(/>) = rContract +-- | Infix synonym for 'scalarProd' +(.|.) = scalarProd +-- | Infix synonym for 'dot' +(<.>) = dot +-- | Infix synonym for 'wedge' +(/\) = wedge + +-- | Left contraction +lContract :: (Eq a, Basis b a) => Clif b a -> Clif b a -> Clif b a +lContract = contractWith (flip (-)) + +-- | Right contraction +rContract :: (Eq a, Basis b a) => Clif b a -> Clif b a -> Clif b a +rContract = contractWith (-) + +-- | Scalar product (0-grade components of the blade products) +scalarProd :: (Eq a, Basis b a) => Clif b a -> Clif b a -> Clif b a +scalarProd = contractWith (const (const 0)) + +-- | Dot product +dot :: (Eq a, Basis b a) => Clif b a -> Clif b a -> Clif b a +dot = contractWith (abs .: (-)) + +-- | Hestenes dot product +hestenes :: (Eq a, Basis b a) => Clif b a -> Clif b a -> Clif b a +hestenes = contractWith f + where f 0 _ = -1 -- negative grades don't exist, so the product vanishes + f _ 0 = -1 + f a b = abs (a - b) + +-- | Wedge product +wedge :: (Eq a, Basis b a) => Clif b a -> Clif b a -> Clif b a +wedge = contractWith (+) + +-- | Hodge dual of a 'Clif' in a Clifford algebra specified by a given pseudoscalar (volume element): +-- +-- prop> hodge (E <$> "abc") $ blade [E 'b'] 1 == blade (E <$> "ac") 1 +-- +hodge :: (Eq a, Basis b a) => [b] -> Clif b a -> Clif b a +hodge = flip (*) . flip blade (-1) + +-- | Projection of Clif x in the direction of Clif y, defined as +-- +-- prop>proj x y == (x <\ recip y) <\ y +-- +proj :: (Eq a, Basis b a, Fractional a) => Clif b a -> Clif b a -> Clif b a +proj x y = (x `lContract` recip y) * y diff --git a/src/Clif/Arbitrary.hs b/src/Clif/Arbitrary.hs new file mode 100644 --- /dev/null +++ b/src/Clif/Arbitrary.hs @@ -0,0 +1,60 @@ +{-| +Module : Clif.Arbitrary +Copyright : (c) Matti A. Eskelinen, 2016-2017 +License : MIT +Maintainer : matti.a.eskelinen@gmail.com +Stability : experimental +Portability : POSIX + +This module provides (orphan) Arbitrary instances and various other generators for creating random 'Clif's and 'Basis' elements using "Test.QuickCheck". + +-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +module Clif.Arbitrary + ( + -- * Generators + ascList, ascLists + , kBlade + ) where +import Clif.Basis +import Clif.Internal + +import Data.List (nub) +import Test.QuickCheck + +-- | Arbitrary instance for 'Euclidean' +instance Arbitrary a => Arbitrary (Euclidean a) where + arbitrary = E <$> arbitrary + shrink = genericShrink + +-- | Arbitrary Instance for 'Lorentzian' +instance Arbitrary a => Arbitrary (Lorentzian a) where + arbitrary = elements [T, S] <*> arbitrary + shrink = genericShrink + +-- | Arbitrary instance for a 'Clif' +instance (Ord b, Arbitrary a, Arbitrary b) => Arbitrary (Clif b a) where + arbitrary = Clif <$> arbitrary + shrink = map Clif . shrink . unClif + +-- | 'orderedList' with only unique elements. +-- Useful for generating blades with e.g. +-- +-- @ +-- 'blade' '<$>' 'ascList' +-- @ +-- +ascList :: (Ord a, Arbitrary a) => Gen [a] +ascList = nub <$> orderedList + +-- | An ascending list split into two at a random point. Useful for generating a pair of blades without common vectors. +ascLists :: (Ord a, Arbitrary a) => Gen ([a], [a]) +ascLists = do + xs <- ascList + split <- elements [0..length xs] + return $ splitAt split xs + +-- | Given k, returns a generator for k-blades ('Clif's containing only a single blade of grade k). +kBlade :: (Eq a, Basis b a, Arbitrary a, Arbitrary b) => Int -> Gen (Clif b a) +kBlade k = blade <$> vector k <*> arbitrary + diff --git a/src/Clif/Basis.hs b/src/Clif/Basis.hs new file mode 100644 --- /dev/null +++ b/src/Clif/Basis.hs @@ -0,0 +1,146 @@ +{-| +Module : Clif.Basis +Copyright : (c) Matti A. Eskelinen, 2016-2017 +License : MIT +Maintainer : matti.a.eskelinen@gmail.com +Stability : experimental +Portability : POSIX + += Motivation +A Clifford algebra is generated by a set of linearly independent (but not necessarily orthogonal) basis vectors and a metric (bilinear form) from the vectors to a field (or unital ring). After the basis and bilinear form are fixed, the Clifford algebra is exactly the quotient of the free algebra of the basis vectors by the equality relation generated by the form (see ). + += Implementation +Defining a basis is abstracted using the two-parameter typeclass 'Basis'. A 'Basis' /b a/ for a Clifford algebra is defined by implementing a 'metric' that takes two vectors of type /b/ (for /basis/) and returns a number of type /a/, which is then used by 'basisMul' to implement multiplication of vectors. For computational reasons, any vector type must have some canonical ordering (i.e. be an instance of 'Ord'). The constructors for 'Euclidean' and 'Lorentzian' types can be used to map any such type with a corresponding metric. + +While not strictly necessary for the abstraction, functions 'canonical' and 'basisMul' are included in the class with respective default implementations ('orderOrthoBasis' and 'orthoMul') in order to provide an easy way to override them using a more efficient implementation for a user-supplied basis. Note that the default implementations are only valid for a diagonal metric; Use 'orderBasis' and 'nonOrthoMul' if you define a non-diagonal metric. + +-} +{-# LANGUAGE + MultiParamTypeClasses + , FlexibleInstances + , DeriveGeneric + #-} +module Clif.Basis + ( -- * Classes + Basis(..) + -- * Basis construction + , Euclidean(..) + , Lorentzian(..) + -- * Multiplication of blades + , orthoMul + , nonOrthoMul + , freeMul + -- * Computation of canonical order + , orderOrthoBasis + , orderBasis + ) where +import GHC.Generics + + +-- | A typeclass for specifying a metric space +-- The default implementation uses 'orthoMul'; Replace it if you require a non-diagonal metric. +class (Ord b, Num a) => Basis b a where + metric :: b -> b -> a + canonical :: ([b], a) -> [([b], a)] + canonical (b, a) = [orderOrthoBasis [] b a] + basisMul :: ([b], a) -> ([b], a) -> [([b], a)] + basisMul a b = [orthoMul a b] + +-- | Data type for constructing Euclidean basis vectors (metric is 1 for matching vectors, 0 otherwise). +newtype Euclidean a = E {getE :: a} + deriving (Eq, Ord, Generic) + +instance Show a => Show (Euclidean a) where + show (E a) = "E " ++ show a + +instance (Ord b, Num a) => Basis (Euclidean b) a where + metric a b = if a == b then 1 else 0 + +-- | Data type for Lorentzian (mixed signature) basis vectors. +data Lorentzian a = T a -- ^ Timelike basis vector (-1) + | S a -- ^ Spacelike basis vector (+1) + deriving (Eq, Ord, Generic) + +instance Show a => Show (Lorentzian a) where + show (T a) = "T " ++ show a + show (S a) = "S " ++ show a + +instance (Ord b, Num a) => Basis (Lorentzian b) a where + metric (T a) (T b) = if a == b then -1 else 0 + metric (S a) (S b) = if a == b then 1 else 0 + metric _ _ = 0 + +-- | Free multiplication of scaled blades (concatenation of vectors and multiplication of scalars). +freeMul :: Num a => ([b], a) -> ([b], a) -> ([b], a) +freeMul (s, m) (t, n) = (s++t, m * n) + +-- | Multiplies together two blades of orthogonal vectors and simplifies the results to canonical order using 'orderOrthoBasis'. +orthoMul :: Basis b a => ([b], a) -> ([b], a) -> ([b], a) +orthoMul = (uncurry (orderOrthoBasis []) .) . freeMul + +-- | Multiplies two blades of (possibly non-orthogonal) basis vectors and simplifies the result to canonical order using 'orderBasis'. +-- Note that this generally results in a direct sum of blades, i.e. a list. +nonOrthoMul :: (Eq a, Basis b a) => ([b], a) -> ([b], a) -> [([b], a)] +nonOrthoMul = (uncurry (orderBasis []) .) . freeMul + +-- | Canonical (ordered) form of a scaled basis blade of orthogonal vectors +-- Uses a gnome sort to commute the basis vectors to order and keeps track of the commutations +-- w.r.t. the blade multiplier, applying the equation +-- +-- /ba = -ab/ +-- +-- for each consecutive pair /ba/ in wrong order. +-- Note that this relation holds only for a basis with a diagonal bilinear form, for others use the more general (albeit slower) 'orderBasis'. +orderOrthoBasis, backOrderOrtho :: Basis b a => [b] -> [b] -> a -> ([b], a) +orderOrthoBasis pre (a:b:rest) c + + -- Vectors are in order: move on. + | a < b = orderOrthoBasis (pre ++ [a]) (b:rest) c + + -- Same vector: simplify using the metric signature for the vector + | a == b = backOrderOrtho pre rest $ metric a b * c + + -- Wrong order: flip and back up one step. Since the vectors are + -- orthogonal, we flip the sign of the geometric product. + | a > b = backOrderOrtho pre (b:a:rest) (-c) +orderOrthoBasis pre rest c = (pre ++ rest, c) + +-- Backwards gnome sort for orthogonal basis +backOrderOrtho [] rest c = orderOrthoBasis [] rest c +backOrderOrtho pre rest c = orderOrthoBasis (init pre) (last pre:rest) c + + +-- | Canonical (ordered) form of a scaled basis blade of vectors +-- Uses a gnome sort to commute the basis vectors to order and keeps track of the commutations +-- applying the equation +-- +-- /ab = 2 B(a,b) - ba/ +-- +-- for any given bilinear form B between the basis vectors and the number ring. +-- Note that this is slower than using 'orderOrthoBasis' if your bilinear form is diagonal (/B(a,b) = 0/ for /a != b/) +orderBasis, backOrder :: (Eq a, Basis b a) => [([b],a)] -> [b] -> a -> [([b],a)] +orderBasis _ _ 0 = [] -- Optimize zero to empty list +orderBasis pre (a:b:rest) c = case compare a b of + -- When in order already, add the first element to each sorted blade + LT -> let add = ([a],1) in orderBasis (if null pre then [add] else map (`freeMul` add) pre) (b:rest) c + -- Simplify equals away and propagate the sort back + EQ -> backOrder pre rest $ metric a b * c + -- When in wrong order, apply the anticommutator relation and branch out + GT -> backOrder pre (b:a:rest) (-c) ++ backOrder pre rest ((metric a b + metric a b) * c) +orderBasis [] rest m = [(rest, m)] +orderBasis pre rest m = map f pre + where f (y, n) = (y ++ rest, n * m) + +-- Gnome sort backwards for non-orthogonal basis +backOrder _ _ 0 = [] +backOrder [] rest c = orderBasis [] rest c +backOrder pre rest c = do + ((prev, s), xs) <- selections pre + case null prev of + True -> orderBasis xs rest (s * c) + False -> orderBasis ((init prev, s):xs) (last prev:rest) c + +-- Given a list, returns all possible tuples containing an element of the original list and the rest of the list. +selections :: [a] -> [(a, [a])] +selections [] = [] +selections (x:xs) = (x,xs):[(y,x:ys) | (y,ys)<-selections xs] diff --git a/src/Clif/Internal.hs b/src/Clif/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Clif/Internal.hs @@ -0,0 +1,190 @@ +{-| +Module : Clif.Internals +Copyright : (c) Matti A. Eskelinen, 2016 +License : MIT +Maintainer : matti.a.eskelinen@gmail.com +Stability : experimental +Portability : POSIX + + +This module implements the various 'Clif' operations on the underlying type. Currently 'Clif' is implemented on top of "Data.Map", with blades as keys and scalar multipliers as values. + += Warning +__This module is not intended to be imported by end users and may change drastically in the future. It is currently exposed (and documented) only to help development.__ + +As development continues, some of the definitions here may be exported from the other modules. Comments and suggestions are welcomed. + +-} +{-# LANGUAGE + GeneralizedNewtypeDeriving + #-} +module Clif.Internal where + +import Clif.Basis + +import Control.Applicative +import qualified Data.Map as M +import Data.List (nub, intercalate) +import Data.Map (Map) +import Data.Function (on) + +-- * Type @Clif@ + +-- | A data type representing a Clif (multivector) composed of direct sum of scaled blades +newtype Clif b a = Clif {unClif :: Map [b] a} + deriving (Functor) + +-- * Constructors + +-- | Constructs a 'Clif' from a list of blades and their multipliers in canonical form. +-- +-- >>> fromList [([], 42), ([E 1, E 2], 1)] +-- 42 *: [] + 1 *: [E 1,E 2] +-- +fromList :: (Eq a, Basis b a) => [([b], a)] -> Clif b a +fromList = Clif . canon' . M.fromListWith (+) + +-- | Constructor for a blade +blade :: (Eq a, Basis b a) => [b] -> a -> Clif b a +blade = Clif .: M.singleton + +-- | Infix synonym for @'flip' 'blade'@ +infix 9 *: +(*:) :: (Eq a, Basis b a) => a -> [b] -> Clif b a +(*:) = flip blade + +-- | Constructor for basis vector values. +-- Note that @'vec' a s@ is equivalent to @'blade' [a] s@. +vec :: (Eq a, Basis b a) => b -> a -> Clif b a +vec = Clif .: M.singleton . (:[]) + +-- | 'Show' instance just shows the underlying 'Map' for now +instance (Show b, Show a) => Show (Clif b a) where + show = intercalate " + " . map showBlade . M.assocs . unClif + where showBlade (b, s) = show s ++ " *: " ++ show b + +-- | The Eq instance calculates the canonical forms of the compared Clifs before comparison. +instance (Eq b, Eq a, Basis b a) => Eq (Clif b a) where + (==) = (==) `on` (unClif . canon) + +-- | Note that abs and signum are only well-defined on the scalar component of each Clif, and zero otherwise. +instance (Eq a, Basis b a) => Num (Clif b a) where + (+) = gPlus + (*) = gMul + abs = fmap abs . grade 0 + negate = fmap negate + signum = fmap signum . grade 0 + fromInteger = blade mempty . fromInteger + +-- | Inverse elements only exist for Clifs c for which c times c is scalar. For others, recip does not terminate. +instance (Eq a, Basis b a, Fractional a) => Fractional (Clif b a) where + fromRational = blade mempty . fromRational + recip v | isScalar v = fmap recip v + | otherwise = revM . recip $ revM v + where revM = gMul (rev v) + +-- * Operations + +-- | The Clifford (geometric) product on 'Clif's. +gMul :: (Eq a, Basis b a) => Clif b a -> Clif b a -> Clif b a +gMul = Clif .: gMul' `on` unClif + +-- | The Clifford product on 'Map's of blades and multipliers. Filter out zero values +gMul' :: (Eq a, Basis b a) => Map [b] a -> Map [b] a -> Map [b] a +gMul' = M.filter (/= 0) . M.fromListWith (+) . concat .: liftA2 basisMul `on` M.assocs + +-- | Addition of 'Clif' values (direct sum). +gPlus :: (Eq a, Basis b a) => Clif b a -> Clif b a -> Clif b a +gPlus = Clif .: gPlus' `on` unClif + +-- | Direct sum of matching keys from two 'Map's. Filter out zero values. +gPlus' :: (Eq a, Basis b a) => Map [b] a -> Map [b] a -> Map [b] a +gPlus' = M.filter (/= 0) .: M.unionWith (+) + +-- | Reverse of a 'Clif', i.e. the reverse of all its component blades. +rev :: Ord b => Clif b a -> Clif b a +rev = Clif . rev' . unClif + +-- | Reverses each blade (key) +rev' :: Ord b => Map [b] a -> Map [b] a +rev' = M.mapKeys reverse + +-- | Returns the canonical form of a 'Clif' +canon :: (Eq a, Basis b a) => Clif b a -> Clif b a +canon = Clif . canon' . unClif + +-- | Returns the canonical representation of a 'Clif' (blades simplified and in canonical order) +canon' :: (Eq a, Basis b a) => Map [b] a -> Map [b] a +canon' = M.filter (/=0) . M.fromListWith (+) . concatMap canonical . M.assocs + +-- | Grade projection on the given grade. For negative values, returns zero. +-- +-- Note that this always calculates the canonical form of a Clif before projecting it. +grade :: (Eq a, Basis b a) => Int -> Clif b a -> Clif b a +grade k = Clif . grade' k . unClif + +-- | Filter blades (keys) by their length. +grade' :: (Eq a, Basis b a) => Int -> Map [b] a -> Map [b] a +grade' k = M.filterWithKey f . canon' + where f = const . (==) k . length + +-- | General product or contraction of Clifs using a given grade function. +-- Given a function @f :: 'Int' -> 'Int' -> 'Int'@ and p, q-grade Clifs A and B, contractWith returns +-- the f(p,q)-grade projection of the 'Clif' A times B. +contractWith :: (Eq a, Basis b a) => (Int -> Int -> Int) -> Clif b a -> Clif b a -> Clif b a +contractWith f = sum .: liftA2 g `on` grades + where g (r, x) (s, y) = let k = f r s in + case compare k 0 of + LT -> 0 + _ -> grade k $ x * y + +-- * Properties + +-- | List of nonzero grades. +-- +-- Note that this always calculates the canonical form of a Clif before recovering the filled grades. +filledGrades :: (Eq a, Basis b a) => Clif b a -> [Int] +filledGrades = filledGrades' . unClif + +-- | List of nonzero grades +-- +-- Note that this always calculates the canonical form of a Clif before recovering the filled grades. +filledGrades' :: (Eq a, Basis b a) => Map [b] a -> [Int] +filledGrades' = nub . map length . M.keys . canon' + +-- | Returns a list containing each non-zero grade component of a 'Clif' and it's grade as an 'Int'. +-- +-- Note that this always calculates the canonical form of a Clif before testing any operations. +grades :: (Eq a, Basis b a) => Clif b a -> [(Int, Clif b a)] +grades v = [(k, grade k v) | k <- filledGrades v] + +-- | True if the 'Clif' contains no nonzero blades of grade greater than zero. +-- +-- Note that this always calculates the canonical form of a Clif before testing whether it is scalar. +isScalar :: (Eq a, Basis b a) => Clif b a -> Bool +isScalar v = ((0==) . maxGrade) v || isZero v + +-- | True for a zero multivector. +-- +-- Note that this always calculates the canonical form of a Clif before testing whether it is zero. +isZero :: (Eq a, Basis b a) => Clif b a -> Bool +isZero = null . canon' . unClif + +-- | The highest nonempty nonzero grade of a Clif. +-- +-- Note that this always calculates the canonical form of a Clif before recovering the highest grade. +maxGrade :: (Eq a, Basis b a) => Clif b a -> Int +maxGrade v | isZero v = 0 + | otherwise = maximum $ filledGrades v + + +-- * Utility functions + +{-# INLINE (.:) #-} +infixl 8 .: +-- |Composition of unary and binary functions, highly useful since two-parameter constructors are ubiquitous here. Redefined here to skip extra dependencies. +-- +-- prop> (f .: g) x y = f (g x y) +-- +(.:) :: (a -> b) -> (c -> d -> a) -> c -> d -> b +(.:) = (.) . (.) diff --git a/src/Clif/Tutorial.hs b/src/Clif/Tutorial.hs new file mode 100644 --- /dev/null +++ b/src/Clif/Tutorial.hs @@ -0,0 +1,137 @@ +{-# OPTIONS_GHC -fno-warn-unused-imports #-} +{-| +Module : Clif.Tutorial +Copyright : (c) Matti A. Eskelinen, 2016-2017 +License : OtherLicense + +/Clifford algebras/ (or /geometric algebras/) are itself mathematically interesting objects, but also a useful tool for vector algebra. This library attempts to make the Clifford algebraic computations easy and at least somewhat computationally efficient, while keeping the implementation as general as possible. + +Since definitions and terminology vary greatly, here is a (non-rigorous) summary of terms that will be used in this tutorial: + + - A /vector/ is just an element of some set (or type). Note that this set may in principle be infinite. + - A /basis/ is a set of vectors along with a bilinear form which maps a pair of them to some number (see "Clif.Basis"). + - A /blade/ is any finite concatenation (/free product/) of vectors which contains each vector at most once, possibly multiplied by a scalar. + - A /scalar/ is just a number, or a number multiplying the /empty product/. + - A 'Clif' is a collection (direct sum) of the empty product, distinct blades and their multipliers (often called a /multivector/). + - The /Clifford algebra/ is the set (type) of 'Clif's with the Clifford (geometric) product and direct summation. + +cf. + + -} +module Clif.Tutorial + ( + -- * Getting started + -- $intro + + -- ** Defining our algebra + -- $basis + + -- ** Construction of 'Clif' values + -- $construction + + -- * Computation + -- $computation + + -- * Copyright + -- $copyright +) +where +import Clif + +{- $copyright +This tutorial is licensed under a + +-} + +{- $intro + To begin, we just need to import the main module of the library, "Clif". + +> import Clif + + This provides us with + + - Constructors for the type 'Clif' with 'Num', 'Eq' and other instances that implement the Clifford algebra, + - Constructors for 'Euclidean' and 'Lorentzian' basis vectors, + - The typeclass 'Basis' for constructing our own algebras, + - Operations of the Clifford algebra defined in "Clif.Algebra". + + +-} + +{- $basis + The type @'Clif' b a@ joins the type @b@ (for /basis/) and some field (or ring) @'Num' a@ together to form a Clifford algebra. Only 'Clif's with matching types @b@ and @a@ can be multiplied directly. To generate the Clifford product between any types @b@ and @a@, we need to specify the bilinear form between them. This is done by providing an instance of the type class @'Basis' b a@. Let us do that for @'Basis' 'Char' Double@: + +>instance Basis Char Double where +> metric 't' 't' = -1 +> metric a b = if a == b then 1 else 0 + + The minimal complete definition for 'Basis' is the function 'metric', which we have here defined to be a diagonal quadratic form on 'Char'. This is all we need to define the reasonably high-dimensional Clifford algebra __Cl__(1,n)(R) where n is the number of Unicode codepoints represented by 'Char', with the signature (-, +, +, ...) for (\'t\', \'a\', \'b\', ...). + + Few notes: + + - We could have used the provided instance for the newtype 'Lorentzian' to wrap 'Char' with a similar metric: + +>instance (Ord b, Num a) => Basis (Lorentzian b) a where +> metric (T a) (T b) = if a == b then -1 else 0 +> metric (S a) (S b) = if a == b then 1 else 0 +> metric _ _ = 0 + + In that case, @'T' a@ for any @'Char' a@ would have the same signature as the character \'t\' in our definition. However, defining a metric for the plain 'Char' is useful for demonstration, since it provides us with pretty printouts. + + - Note that as in the definition for @'Basis' ('Lorentzian' b) a@, we do not actually need to fix the field @ə@ apart from the 'Num' constraint while defining the basis. + + - If the metric is not diagonal (@'metric' a b /= 0@ for some @a /= b@), we need to replace the default implementations of the functions 'canonical' and 'basisMul' in the instance with more general implementations. See "Clif.Basis" for details. + +-} + + +{- $construction + + Using our 'Char' basis, we can start constructing values of type @'Clif' 'Char' 'Double'@ using the provided constructors. We can start by introducing the vectors needed to represent __Cl__(1,3)(R): + +> t = vec 't' 1 +> x = vec 'x' 1 +> y = vec 'y' 1 +> z = vec 'z' 1 + + To make it specific that we are working in __Cl__(1,3)(R), we can define the /pseudoscalar/ @txyz@. All the following definitions are equivalent: + +> i = t * x * y * z + +> i = blade "txyz" 1 + +> i = 1 *: "txyz" + + Note that the last form using the infix operator '(*:)' is used for the 'Show' instance to produce concise output. + p + Due to the 'Num' instance, we do not usually need to explicitly embed scalars. If we want to be specific, we can write + +> answer = 42 :: Clif Char Double + +-} + +{- $computation + + We can now use any of the available operations to calculate on the 'Clif' values, such as + + - Simple multivector algebra: + +>>> 2 * x * y - y * x +1.0 *: "xy" + + - Wedge products: + +>>> x /\ y +1.0 *: "xy" + + - Reversion: + +>>> rev i +1.0 *: "zyxt" + + - Since 'Double's have inverses, so do 'Clif's for which @v * v@ is scalar. Trivial example is vectors: + +>>> x * y * z / y +-1.0 *: "xy" + +-} diff --git a/tests/tests.hs b/tests/tests.hs new file mode 100644 --- /dev/null +++ b/tests/tests.hs @@ -0,0 +1,19 @@ +module Main where +import Test.Tasty +import qualified BasisTests as BT +import qualified ClifTests as CT +import qualified InternalTests as IT + +main :: IO () +main = defaultMain tests + +tests :: TestTree +tests = testGroup "Tests" + [ testGroup "Implementation" [IT.tests] + , testGroup "Properties" + [testGroup "Basis" [BT.tests] + ,testGroup "Clifford algebra" [CT.tests] + ] + ] + +