@@ -0,0 +1,254 @@
1+module Data.Matrix
2+
3+import Data.Complex
4+import Data.ZZ
5+
6+%default total
7+
8+infixr 2 <.> -- vector inner product
9+infixr 2 >< -- vector outer product
10+infixr 2 <<>> -- matrix commutator
11+infixr 2 >><< -- matrix anticommutator
12+infixl 3 <\> -- row times a matrix
13+infixr 4 </> -- matrix times a column
14+infixr 5 <> -- matrix multiplication
15+infixr 7 \&\ -- vector tensor product
16+infixr 7 <&> -- matrix tensor product
17+
18+-----------------------------------------------------------------------
19+-- Vectors as members of algebraic classes
20+-----------------------------------------------------------------------
21+
22+instance Semigroup a => Semigroup (Vect n a) where
23+ (<+>) v w = zipWith (<+>) v w
24+
25+instance Monoid a => Monoid (Vect n a) where
26+ neutral {n} = replicate n neutral
27+
28+instance Group a => Group (Vect n a) where
29+ inverse = map inverse
30+
31+instance AbelianGroup a => AbelianGroup (Vect n a) where {}
32+
33+instance Ring a => Ring (Vect n a) where
34+ (<*>) v w = zipWith (<*>) v w
35+
36+instance RingWithUnity a => RingWithUnity (Vect n a) where
37+ unity {n} = replicate n unity
38+
39+instance Field a => Field (Vect n a) where
40+ inverseM = map inverseM
41+
42+instance RingWithUnity a => Module a (Vect n a) where
43+ (<#>) r v = map (r <*>) v
44+
45+instance RingWithUnity a => Module a (Vect n (Vect l a)) where
46+ (<#>) r m = map (r <#>) m
47+-- should be Module a b => Module a (Vect n b), but results in 'overlapping instance'
48+
49+-----------------------------------------------------------------------
50+-- (Ring) Vector functions
51+-----------------------------------------------------------------------
52+
53+||| Inner product of ring vectors
54+(<.>) : Ring a => Vect n a -> Vect n a -> a
55+(<.>) w v = foldr (<+>) neutral (zipWith (<*>) w v)
56+
57+||| Tensor multiply (⊗) ring vectors
58+(\&\) : Ring a => Vect n a -> Vect m a -> Vect (n * m) a
59+(\&\) {n} {m} v w = zipWith (<*>) (oextend m v) (orep n w) where
60+ orep : (n : Nat) -> Vect m a -> Vect (n * m) a
61+ orep n v = concat $ replicate n v
62+ oextend : (n : Nat) -> Vect k a -> Vect (k * n) a
63+ oextend n w = concat $ map (replicate n) w
64+
65+||| Standard basis vector with one nonzero entry, ring data type and vector-length unfixed
66+basis : RingWithUnity a => {d : Nat} -> (Fin d) -> Vect d a
67+basis i = replaceAt i unity $ neutral
68+
69+-----------------------------------------------------------------------
70+-- Matrix functions
71+-----------------------------------------------------------------------
72+
73+||| Matrix with n rows and m columns
74+Matrix : Nat -> Nat -> Type -> Type
75+Matrix n m a = Vect n (Vect m a)
76+
77+||| Gets the specified column of a matrix. For rows use the vector function 'index'
78+getCol : Fin m -> Matrix n m a -> Vect n a
79+getCol fm q = map (index fm) q
80+
81+||| Deletes the specified column of a matrix. For rows use the vector function 'deleteAt'
82+deleteCol : Fin (S m) -> Matrix n (S m) a -> Matrix n m a
83+deleteCol f m = map (deleteAt f) m
84+
85+||| Matrix element at specified row and column indices
86+indices : Fin n -> Fin m -> Matrix n m a -> a
87+indices f1 f2 m = index f2 (index f1 m)
88+
89+||| Matrix times a column vector
90+(</>) : Ring a => Matrix n m a -> Vect m a -> Vect n a
91+(</>) m v = map (v <.>) m
92+
93+||| Matrix times a row vector
94+(<\>) : Ring a => Vect n a -> Matrix n m a -> Vect m a
95+(<\>) r m = map (r <.>) $ transpose m
96+
97+||| Matrix multiplication
98+(<>) : Ring a => Matrix n k a ->
99+ Matrix k m a ->
100+ Matrix n m a
101+(<>) m1 m2 = map (<\> m2) m1
102+
103+||| Tensor multiply (⊗) for ring matrices
104+(<&>) : Ring a => Matrix h1 w1 a -> Matrix h2 w2 a -> Matrix (h1 * h2) (w1 * w2) a
105+(<&>) m1 m2 = zipWith (\&\) (stepOne m1 m2) (stepTwo m1 m2) where
106+ stepOne : Matrix h1 w1 a -> Matrix h2 w2 a -> Matrix (h1 * h2) w1 a
107+ stepOne {h2} m1 m2 = concat $ map (replicate h2) m1
108+ stepTwo : Matrix h1 w1 a -> Matrix h2 w2 a -> Matrix (h1 * h2) w2 a
109+ stepTwo {h1} m1 m2 = concat $ (Vect.replicate h1) m2
110+
111+||| Cast a vector from a standard Vect to a proper n x 1 matrix
112+col : Vect n a -> Matrix n 1 a
113+col v = map (\x => [x]) v
114+
115+||| Cast a row from a standard Vect to a proper 1 x n matrix
116+row : Vect n a -> Matrix 1 n a
117+row r = [r]
118+
119+||| Outer product between ring vectors
120+(><) : Ring a => Vect n a -> Vect m a -> Matrix n m a
121+(><) x y = (col x) <> (row y)
122+
123+||| All finite numbers up to the given bound
124+allN : (n : Nat) -> Vect n (Fin n)
125+allN Z = Nil
126+allN (S n) = FZ :: (map FS $ allN n)
127+
128+||| Identity matrix
129+Id : RingWithUnity a => {d : Nat} -> Matrix d d a
130+Id {d} = map (\n => basis n) $ allN d
131+
132+||| Matrix commutator
133+(<<>>) : Ring a => Matrix n n a -> Matrix n n a -> Matrix n n a
134+(<<>>) m1 m2 = (m1 <> m2) <-> (m2 <> m1)
135+
136+||| Matrix anti-commutator
137+(>><<) : Ring a => Matrix n n a -> Matrix n n a -> Matrix n n a
138+(>><<) m1 m2 = (m1 <> m2) <+> (m2 <> m1)
139+
140+-----------------------------------------------------------------------
141+-- Matrix Algebra Properties
142+-----------------------------------------------------------------------
143+
144+-- TODO: Prove properties of matrix algebra for 'Verified' algebraic classes
145+
146+-----------------------------------------------------------------------
147+-- Numberic data types as rings
148+-----------------------------------------------------------------------
149+
150+instance Semigroup Integer where
151+ (<+>) = (+)
152+
153+instance Monoid Integer where
154+ neutral = 0
155+
156+instance Group Integer where
157+ inverse = (* -1)
158+
159+instance AbelianGroup Integer
160+
161+instance Ring Integer where
162+ (<*>) = (*)
163+
164+instance RingWithUnity Integer where
165+ unity = 1
166+
167+
168+instance Semigroup Int where
169+ (<+>) = (+)
170+
171+instance Monoid Int where
172+ neutral = 0
173+
174+instance Group Int where
175+ inverse = (* -1)
176+
177+instance AbelianGroup Int
178+
179+instance Ring Int where
180+ (<*>) = (*)
181+
182+instance RingWithUnity Int where
183+ unity = 1
184+
185+
186+instance Semigroup Float where
187+ (<+>) = (+)
188+
189+instance Monoid Float where
190+ neutral = 0
191+
192+instance Group Float where
193+ inverse = (* -1)
194+
195+instance AbelianGroup Float
196+
197+instance Ring Float where
198+ (<*>) = (*)
199+
200+instance RingWithUnity Float where
201+ unity = 1
202+
203+instance Field Float where
204+ inverseM f = 1 / f
205+
206+
207+instance Semigroup Nat where
208+ (<+>) = (+)
209+
210+instance Monoid Nat where
211+ neutral = 0
212+
213+instance VerifiedSemigroup Nat where
214+ semigroupOpIsAssociative = plusAssociative
215+
216+instance VerifiedMonoid Nat where
217+ monoidNeutralIsNeutralL = plusZeroRightNeutral
218+ monoidNeutralIsNeutralR = plusZeroLeftNeutral
219+
220+
221+instance Semigroup ZZ where
222+ (<+>) = (+)
223+
224+instance Monoid ZZ where
225+ neutral = 0
226+
227+instance Group ZZ where
228+ inverse = (* -1)
229+
230+instance AbelianGroup ZZ
231+
232+instance Ring ZZ where
233+ (<*>) = (*)
234+
235+instance RingWithUnity ZZ where
236+ unity = 1
237+
238+
239+instance Semigroup a => Semigroup (Complex a) where
240+ (<+>) (a :+ b) (c :+ d) = (a <+> c) :+ (b <+> d)
241+
242+instance Monoid a => Monoid (Complex a) where
243+ neutral = (neutral :+ neutral)
244+
245+instance Group a => Group (Complex a) where
246+ inverse (r :+ i) = (inverse r :+ inverse i)
247+
248+instance Ring a => AbelianGroup (Complex a) where {}
249+
250+instance Ring a => Ring (Complex a) where
251+ (<*>) (a :+ b) (c :+ d) = (a <*> c <-> b <*> d) :+ (a <*> d <+> b <*> c)
252+
253+instance RingWithUnity a => RingWithUnity (Complex a) where
254+ unity = (unity :+ neutral)