symgeoalg

symbolic geometric algebra calculations
git clone git://git.rr3.xyz/symgeoalg
Log | Files | Refs | README | LICENSE

commit ab38267ed2c9b057de5be6780c0e9184be28dd0b
parent 884f33c2c9174935dfa5e258bdba4accb40c7d10
Author: Robert Russell <robertrussell.72001@gmail.com>
Date:   Fri, 26 Jul 2024 23:52:18 -0700

Organize and save

Currently reworking GeometricAlgebra.

I just want to get some stuff in git history before I delete it.

Diffstat:
MMain.lean | 410++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------
1 file changed, 283 insertions(+), 127 deletions(-)

diff --git a/Main.lean b/Main.lean @@ -1,5 +1,45 @@ import Batteries.Data.Rat +-------------------------------------------------------------------------------- +-- Util + +def finSum (i : Fin (m + n)) : Sum (Fin m) (Fin n) := + -- Commute m and n in i's type to please Fin.subNat. + have fin_comm : Fin (m + n) = Fin (n + m) := by rw [Nat.add_comm] + let i : Fin (n + m) := cast fin_comm i + + if h : i.val < m then + Sum.inl (Fin.castLT i h) + else + have m_le_i : m <= i.val := Nat.le_of_not_gt h + Sum.inr (Fin.subNat m i m_le_i) + +def finProd (i : Fin (m * n)) : Prod (Fin m) (Fin n) := + let mul_gt_zero_imp_gt_zero {m n : Nat} : (m * n > 0) -> (m > 0) := by + intro + rw [<-Nat.zero_div n] + apply Nat.div_lt_of_lt_mul + rw [Nat.mul_comm] + assumption + + have mn_gt_zero : m * n > 0 := Fin.pos i + have m_gt_zero : m > 0 := mul_gt_zero_imp_gt_zero mn_gt_zero + let q := Fin.ofNat' (i.val / n) m_gt_zero + + have nm_gt_zero : n * m > 0 := by rw [Nat.mul_comm]; exact mn_gt_zero + have n_gt_zero : n > 0 := mul_gt_zero_imp_gt_zero nm_gt_zero + let r := Fin.ofNat' (i.val % n) n_gt_zero + + (q, r) + +-------------------------------------------------------------------------------- +-- Notation + +instance : HPow Type Nat Type where + hPow Y n := Fin n -> Y + +def vec (coords : List Y) : Y ^ coords.length := coords.get + -- SMul represents the ability to use "⋅" as an infix operator. class SMul (α β : Type) where smul : α -> β -> β @@ -8,6 +48,12 @@ class SMul (α β : Type) where -- How do I figure out precedences of standard Lean operators? infixr:75 " ⋅ " => SMul.smul +class Inv (α : Type) where + inv : α -> α + +-------------------------------------------------------------------------------- +-- Single additive operator algebraic structures + class AddSemigroup (A : Type) where add : A -> A -> A class AddMonoid (A : Type) extends AddSemigroup A where @@ -16,7 +62,7 @@ class AddGroup (A : Type) extends AddMonoid A where sub : A -> A -> A neg : A -> A := fun x => sub zero x --- Notation for AddXXX +-- Notation type class instances -- We could just, e.g., replace AddSemigroup A with Add A, but I think we ought -- to distinguish "notation type classes" from "conceptual type classes". It's -- especially important in this case, because these structures should satisfy @@ -33,6 +79,27 @@ instance [AddGroup A] : Sub A where instance [AddGroup A] : Neg A where neg := AddGroup.neg +instance [AddSemigroup A] [AddSemigroup B] : AddSemigroup (A × B) where + add x y := (x.1 + y.1, x.2 + y.2) +instance [AddMonoid A] [AddMonoid B] : AddMonoid (A × B) where + zero := (0, 0) +instance [AddGroup A] [AddGroup B] : AddGroup (A × B) where + sub x y := (x.1 - y.1, x.2 - y.2) + +instance AddSemigroup.instFun {X : Type} [AddSemigroup A] : AddSemigroup (X -> A) where + add f g := fun x => f x + g x +instance AddMonoid.instFun {X : Type} [AddMonoid A] : AddMonoid (X -> A) where + zero := fun _ => 0 +instance AddGroup.instFun {X : Type} [AddGroup A] : AddGroup (X -> A) where + sub f g := fun x => f x - g x + +instance [AddSemigroup A] {n : Nat} : AddSemigroup (A ^ n) := AddSemigroup.instFun +instance [AddMonoid A] {n : Nat} : AddMonoid (A ^ n) := AddMonoid.instFun +instance [AddGroup A] {n : Nat} : AddGroup (A ^ n) := AddGroup.instFun + +-------------------------------------------------------------------------------- +-- Single multiplicative operator algebraic structures + class MulSemigroup (M : Type) where mul : M -> M -> M class MulMonoid (M : Type) extends MulSemigroup M where @@ -41,36 +108,71 @@ class MulGroup (M : Type) extends MulMonoid M where div : M -> M -> M inv : M -> M := fun x => div one x --- Notation for MulXXX +-- Notation type class instances instance [MulMonoid M] : OfNat M 1 where ofNat := MulMonoid.one instance [MulSemigroup M] : Mul M where mul := MulSemigroup.mul instance [MulGroup M] : Div M where div := MulGroup.div -def inv := @MulGroup.inv +instance [MulGroup M] : Inv M where + inv := MulGroup.inv + +instance [MulSemigroup A] [MulSemigroup B] : MulSemigroup (A × B) where + mul x y := (x.1 * y.1, x.2 * y.2) +instance [MulMonoid A] [MulMonoid B] : MulMonoid (A × B) where + one := (1, 1) +instance [MulGroup A] [MulGroup B] : MulGroup (A × B) where + div x y := (x.1 / y.1, x.2 / y.2) + +instance MulSemigroup.instFun {X : Type} [MulSemigroup M] : MulSemigroup (X -> M) where + mul f g := fun x => f x * g x +instance MulMonoid.instFun {X : Type} [MulMonoid M] : MulMonoid (X -> M) where + one := fun _ => 1 +instance MulGroup.instFun {X : Type} [MulGroup M] : MulGroup (X -> M) where + div f g := fun x => f x / g x + +instance [MulSemigroup M] {n : Nat} : MulSemigroup (M ^ n) := MulSemigroup.instFun +instance [MulMonoid M] {n : Nat} : MulMonoid (M ^ n) := MulMonoid.instFun +instance [MulGroup M] {n : Nat} : MulGroup (M ^ n) := MulGroup.instFun + +-------------------------------------------------------------------------------- +-- Fields +-- TODO: Generalize Field/VectorSpace to Ring/Module? -- Convention: (∀a. div a 0 = 0) and (inv 0 = 0). class Field (F : Type) extends AddGroup F, MulGroup F +instance : Field Rat where + zero := 0 + add := Rat.add + sub := Rat.sub + one := 1 + mul := Rat.mul + div := Rat.div + +-------------------------------------------------------------------------------- +-- Vector spaces + -- F is an outParam, because usually V will only be a vector space over one -- field, and this allows us to leave F implicit in many cases. class VectorSpace (F : outParam Type) [Field F] (V : Type) [AddGroup V] where smul : F -> V -> V --- Notation for VectorSpace instance [Field F] [AddGroup V] [VectorSpace F V] : SMul F V where smul := VectorSpace.smul --- Fields F are vector spaces over F. instance [Field F] : VectorSpace F F where smul := MulSemigroup.mul --- We could define Coordinates V n = NSpace n F, but I think it's important to --- recognize the conceptual difference between these two types: Coordinates V n --- is necessarily a tuple of scalars, with no further inherent structure, --- whereas NSpace n F is a finite product of vector spaces, which just happens --- to look like a tuple of scalars. +instance [Field F] [AddGroup V] [VectorSpace F V] [AddGroup W] [VectorSpace F W] : VectorSpace F (V × W) where + smul a p := (a ⋅ p.1, a ⋅ p.2) + +instance VectorSpace.instFun [Field F] [AddGroup V] [VectorSpace F V] : VectorSpace F (X -> V) where + smul a f := fun x => a ⋅ f x + +instance [Field F] [AddGroup V] [VectorSpace F V] {n : Nat} : VectorSpace F (V ^ n) := VectorSpace.instFun + def Coordinates [Field F] (V : Type) [AddGroup V] [VectorSpace F V] (n : Nat) : Type := Fin n -> F @@ -90,11 +192,33 @@ def FiniteBasis.vector [Field F] [AddGroup V] [VectorSpace F V] (b : FiniteBasis V n) (coords : Coordinates V n) : V := Fin.foldl n (fun acc i => acc + coords i ⋅ b.basis i) 0 -def fieldBasis [Field F] : FiniteBasis F 1 := +def FiniteBasis.field (F : Type) [Field F] : FiniteBasis F 1 := ⟨fun 0 => 1, fun a 0 => a⟩ --- OrthoQuadraticForm V B is a quadratic form on V with respect to which the --- finite basis B is orthogonal. +def FiniteBasis.mul [Field F] [AddGroup V] [VectorSpace F V] (b : FiniteBasis V m) + [AddGroup W] [VectorSpace F W] (c : FiniteBasis W n) : FiniteBasis (V × W) (m + n) := + let basis i := + match finSum i with + | .inl i => (b.basis i, 0) + | .inr i => (0, c.basis i) + let coords v i := + match finSum i with + | .inl i => b.coords v.1 i + | .inr i => c.coords v.2 i + ⟨basis, coords⟩ + +def FiniteBasis.pow [Field F] [AddGroup V] [VectorSpace F V] (b : FiniteBasis V m) + (n : Nat) : FiniteBasis (V ^ n) (m * n) := + let basis i := + let ⟨iq, ir⟩ := finProd i + fun j => if j = ir then b.basis iq else 0 + let coords v i := + let ⟨iq, ir⟩ := finProd i + b.coords (v ir) iq + ⟨basis, coords⟩ + +-- OrthoQuadraticForm b is a quadratic form with respect to which the finite +-- basis b is orthogonal. structure OrthoQuadraticForm [Field F] [AddGroup V] [VectorSpace F V] (_ : FiniteBasis V n) where evalBasis : Fin n -> F @@ -107,146 +231,178 @@ def OrthoQuadraticForm.eval [Field F] [AddGroup V] [VectorSpace F V] -- = ∑ (vCoords i)^2 * q.evalBasis i Fin.foldl n (fun acc i => acc + vCoords i * vCoords i * q.evalBasis i) 0 ---class Algebra (F : Type) [Field F] (A : Type) [AddGroup A] [MulMonoid A] extends VectorSpace F A +-------------------------------------------------------------------------------- +-- Algebras (over a field) + +class Algebra (F : outParam Type) [Field F] (A : Type) [AddGroup A] [MulMonoid A] extends VectorSpace F A + +instance [Field F] : Algebra F F where + +-- TODO: More Algebra instances + + + + +def bitVecFoldl (f : α -> Bool -> α) (bv : BitVec n) (init : α) : α := + sorry + +def bitVecEnuml (n : Nat) (f : α -> BitVec n -> α) (init : α) : α := + let loop acc i := f acc (BitVec.ofFin i) + Fin.foldl (2^n) loop init + --- TODO: Rename to FiniteBasisSet? Usually "Blade" means something more general. -structure Blade (n : Nat) : Type where +structure BasisBlade (n : Nat) : Type where bits : BitVec n -def Blade.has (b : Blade n) (i : Fin n) : Bool := b.bits.getLsb i +-- TODO: ? +def BasisBlade.empty : BasisBlade n -> Bool := sorry + +def BasisBlade.has (bld : BasisBlade n) (i : Fin n) : Bool := + bld.bits.getLsb i + +-- TODO: delete +-- def BasisBlade.cons (bit : Bool) (bld : BasisBlade i) : BasisBlade (i + 1) := +-- ⟨BitVec.cons bit bld.bits⟩ +-- def BasisBlade.enum (n : Nat) : List (BasisBlade n) := +-- let enum (i : Fin n) (b : List (BasisBlade i.val)) : List (BasisBlade (i.val + 1)) := +-- (BasisBlade.cons false <$> b) ++ (BasisBlade.cons true <$> b) +-- Fin.hIterate (List ∘ BasisBlade) [⟨BitVec.nil⟩] enum + +def BasisBlade.mulScalar [Field F] [AddGroup V] [VectorSpace F V] {b : FiniteBasis V n} + (q : OrthoQuadraticForm b) (bld0 bld1 : BasisBlade n) : F := + bitVecFoldl x (bld0 &&& bld1) 0 + +def BasisBlade.foldl (n : Nat) (f : α -> BasisBlade n -> α) (init : α) : α := + let loop acc i := f acc ⟨BitVec.ofFin i⟩ + Fin.foldl (2^n) loop init + +instance : ToString (BasisBlade n) where + toString bld := + let indexToString i acc := if bld.has i then toString i :: acc else acc + let indexStrings := Fin.foldr n indexToString [] + let sep := if n > 9 then "," else "" + "{" ++ String.intercalate sep indexStrings ++ "}" + -def Blade.cons (bit : Bool) (b : Blade i) : Blade (i + 1) := ⟨BitVec.cons bit b.bits⟩ -def Blade.enum (n : Nat) : List (Blade n) := - let enum (i : Fin n) (b : List (Blade i.val)) : List (Blade (i.val + 1)) := - (Blade.cons false <$> b) ++ (Blade.cons true <$> b) - Fin.hIterate (List ∘ Blade) [⟨BitVec.nil⟩] enum -instance : ToString (Blade n) where - toString b := - let vectorIndexStrings := Fin.foldl n (fun acc i => if b.has i then toString (n - i) :: acc else acc) [] - if b.bits != 0 then - "e" ++ String.join vectorIndexStrings - else - "" structure GeometricAlgebra [Field F] [AddGroup V] [VectorSpace F V] {b : FiniteBasis V n} (q : OrthoQuadraticForm b) : Type where - fn : FunctionSpace (Blade n) F + coords : BitVec n -> F + +instance [Field F] [DecidableEq F] [ToString F] [AddGroup V] [VectorSpace F V] + {b : FiniteBasis V n} {q : OrthoQuadraticForm b} : ToString (GeometricAlgebra q) where + toString v := + let bldToString (bld : BitVec n) := + let indexToString (i : Fin n) acc := + if bld.getLsb i then toString i :: acc else acc + let indexStrings := Fin.foldr n indexToString [] + let sep := if n > 9 then "," else "" + "{" ++ String.intercalate sep indexStrings ++ "}" + + let termToString acc bld := + if v.coords bld != 0 then + let coeffString := toString (v.coords bld) + let termString := + if bld != 0 then + s!"{coeffString} * {bldToString bld}" + else + coeffString + termString :: acc + else + acc + let termStrings := bitVecEnuml n termToString [] + + if termStrings.length > 0 then + String.intercalate " + " termStrings + else + "0" instance [Field F] [AddGroup V] [VectorSpace F V] {b : FiniteBasis V n} {q : OrthoQuadraticForm b} : AddGroup (GeometricAlgebra q) where zero := ⟨0⟩ - add u v := ⟨u.fn + v.fn⟩ - sub u v := ⟨u.fn - v.fn⟩ + add v w := ⟨v.coords + w.coords⟩ + sub v w := ⟨v.coords - w.coords⟩ instance [Field F] [AddGroup V] [VectorSpace F V] {b : FiniteBasis V n} {q : OrthoQuadraticForm b} : VectorSpace F (GeometricAlgebra q) where - smul a v := ⟨a ⋅ v.fn⟩ - --- TODO: Rename variables. --- Given u,v : GeometricAlgebra F n and b : Blade n, the b-th component of the --- product u * v is a sum of 2^n terms, where each term is the product of the --- l-th component of u and the r-th component of v, for all l,r : Blade n that multiply to b (modulo a scalar). The type + smul a v := ⟨a ⋅ v.coords⟩ + +-- Multiplication algorithm: +-- Given v w : GeometricAlgebra q and bld : BasisBlade n, the bld-th component +-- of the product v * w is a sum of 2^n terms, where each term is the product of +-- the vbld-th component of v and the wbld-th component of w, for all +-- vbld wbld : BasisBlade n that multiply (XOR) to bld (modulo a scalar). For +-- example, if we take n = 3 and denote blades as bit vectors, then we can write +-- generic vectors v and w as follows: +-- v0 * 000 + v1 * 001 + v2 * 010 + v3 * 011 + v4 * 100 + v5 * 101 + v6 * 110 + v7 * 111 +-- w0 * 000 + w1 * 001 + w2 * 010 + w3 * 011 + w4 * 100 + w5 * 101 + w6 * 110 + w7 * 111 +-- Now, if we want the component of v * w corresponding to blade bld = 010, we +-- consider the following pairs of vbld and wbld: +-- vbld wbld +-- 010 000 +-- 011 001 +-- 000 010 +-- 001 011 +-- 110 100 +-- 111 101 +-- 100 110 +-- 101 111 +-- (These are just all 2^3 3-bit numbers, but with one column (vbld in this +-- case) XORed with bld.) Then we multiply the coefficients corresponding to +-- each row and sum the results. We call each row a "term". We build up each +-- row (i.e., (vbld, wbld) pair) incrementally, bit-by-bit, so an + +-- When vbld and wbld share a bit, it "cancels" +-- and gets replaced by a scalar in accordance with the quadratic form. Also, + + +-- The type -- PartialTerm F i represents part of such a term, specifically the first i bits of l and r -- and an associated scalar a. We also keep track of the parity of the l bits seen so far, so -- that we know if we should flip the sign of a when we get the next r bit. This type should -- be local to the multiplication function, but I don't think structures can be local in Lean -- (ughh!). -structure PartialTerm (F : Type) [Field F] (i : Nat) where - a : F - lParity : F - l : BitVec i - r : BitVec i - -def PartialTerm.cons [Field F] (a : F) (l r : Bool) (pt : PartialTerm F i) : PartialTerm F (i + 1) := - let a' := a * pt.a * if r then pt.lParity else 1 - let lParity' := pt.lParity * if l then -1 else 1 - ⟨a', lParity', BitVec.cons l pt.l, BitVec.cons r pt.r⟩ +-- structure PartialTerm (F : Type) [Field F] (i : Nat) where +-- a : F +-- lParity : F +-- l : BitVec i +-- r : BitVec i +-- def PartialTerm.cons [Field F] (a : F) (l r : Bool) (pt : PartialTerm F i) : PartialTerm F (i + 1) := +-- let a' := a * pt.a * if r then pt.lParity else 1 +-- let lParity' := pt.lParity * if l then -1 else 1 +-- ⟨a', lParity', BitVec.cons l pt.l, BitVec.cons r pt.r⟩ +-- def PartialTerm.eval [Field F] {q : QuadraticForm F n} (u v : GeometricAlgebra F n q) (pt : PartialTerm F n) : F := +-- pt.a * u.fn ⟨pt.l⟩ * v.fn ⟨pt.r⟩ -def PartialTerm.eval [Field F] {q : QuadraticForm F n} (u v : GeometricAlgebra F n q) (pt : PartialTerm F n) : F := - pt.a * u.fn ⟨pt.l⟩ * v.fn ⟨pt.r⟩ - -instance [Field F] {q : QuadraticForm F n} : MulMonoid (GeometricAlgebra F n q) where +instance [Field F] [AddGroup V] [VectorSpace F V] {b : FiniteBasis V n} + {q : OrthoQuadraticForm b} : MulMonoid (GeometricAlgebra q) where one := ⟨fun b => if b.bits = 0 then 1 else 0⟩ - mul u v := - let terms (b : Blade n) (i : Fin n) (acc : List (PartialTerm F i.val)) : List (PartialTerm F (i.val + 1)) := - if b.has i then - (PartialTerm.cons 1 false true <$> acc) ++ (PartialTerm.cons 1 true false <$> acc) - else - (PartialTerm.cons 1 false false <$> acc) ++ (PartialTerm.cons (q i) true true <$> acc) - let fn (b : Blade n) : F := - let ts := Fin.hIterate (List ∘ PartialTerm F) [⟨1, 1, BitVec.nil, BitVec.nil⟩] (terms b) - let x := PartialTerm.eval u v <$> ts - List.foldl AddSemigroup.add 0 x - ⟨fn⟩ - --- TODO: Generalize this printing functionality to arbitrary FunctionSpaces with --- enumerable domain? -instance [Field F] [ToString F] {q : QuadraticForm F n} : ToString (GeometricAlgebra F n q) where - toString v := - let termToString (b : Blade n) : String := - toString (v.fn b) ++ if b.bits != 0 then s!" * {toString b}" else "" - String.intercalate " + " (termToString <$> Blade.enum n) + mul v w := + let coords (bld : BasisBlade n) : F := + let addTerm (acc : F) (vbld : BasisBlade n) : F := + let wbld : BasisBlade n := ⟨vbld.bits ^^^ bld.bits⟩ + let bldScalar : F := + let vCoeff : F := v.coords vbld + let wCoeff : F := w.coords wbld + acc + v.coords vbld * w.coords wbld * vbld.mul wbld + BasisBlade.foldl n addTerm 0 + ⟨coords⟩ + --mul v w := + -- let terms (b : Blade n) (i : Fin n) (acc : List (PartialTerm F i.val)) : List (PartialTerm F (i.val + 1)) := + -- if b.has i then + -- (PartialTerm.cons 1 false true <$> acc) ++ (PartialTerm.cons 1 true false <$> acc) + -- else + -- (PartialTerm.cons 1 false false <$> acc) ++ (PartialTerm.cons (q i) true true <$> acc) + -- let coords (b : BasisBlade n) : F := + -- let ts := Fin.hIterate (List ∘ PartialTerm F) [⟨1, 1, BitVec.nil, BitVec.nil⟩] (terms b) + -- let x := PartialTerm.eval u v <$> ts + -- List.foldl AddSemigroup.add 0 x + -- ⟨coords⟩ -------------------------------------------------------------------------------- --- TODO: Remaining instances for products. -instance [AddSemigroup A] [AddSemigroup B] : AddSemigroup (A × B) where - add x y := (x.1 + y.1, x.2 + y.2) -instance [AddMonoid A] [AddMonoid B] : AddMonoid (A × B) where - zero := (0, 0) -instance [AddGroup A] [AddGroup B] : AddGroup (A × B) where - sub x y := (x.1 - y.1, x.2 - y.2) - -instance : Field Rat where - zero := 0 - add := Rat.add - sub := Rat.sub - one := 1 - mul := Rat.mul - div := Rat.div - --- Given a type X and a vector space V over F, the type of functions X -> V is a --- vector space over F. -def FunctionSpace (X V : Type) : Type := X -> V -instance [AddGroup V] : AddGroup (FunctionSpace X V) where - zero := fun _ => 0 - add f g := fun x => f x + g x - sub f g := fun x => f x - g x -instance [Field F] [AddGroup V] [VectorSpace F V] : VectorSpace F (FunctionSpace X V) where - smul a f := fun x => a ⋅ f x - --- Function spaces with finite domain -abbrev NSpace (n : Nat) (V : Type) : Type := FunctionSpace (Fin n) V - -def divmod (i : Fin (m * n)) : Fin m × Fin n := - let mul_gt_zero_imp_gt_zero {m n : Nat} : (m * n > 0) -> (m > 0) := by - intro - rw [<-Nat.zero_div n] - apply Nat.div_lt_of_lt_mul - rw [Nat.mul_comm] - assumption - have mn_gt_zero : m * n > 0 := Fin.pos i - have m_gt_zero : m > 0 := mul_gt_zero_imp_gt_zero mn_gt_zero - have nm_gt_zero : n * m > 0 := by rw [Nat.mul_comm]; exact mn_gt_zero - have n_gt_zero : n > 0 := mul_gt_zero_imp_gt_zero nm_gt_zero - let q := Fin.ofNat' (i.val / n) m_gt_zero - let r := Fin.ofNat' (i.val % n) n_gt_zero - (q, r) - -def repeatBasis [Field F] [AddGroup V] [VectorSpace F V] (b : FiniteBasis V m) - (n : Nat) : FiniteBasis (NSpace n V) (m * n) := - let basis i := - let ⟨iq, ir⟩ := divmod i - fun j => if j = ir then b.basis iq else 0 - let coords v i := - let ⟨iq, ir⟩ := divmod i - b.coords (v ir) iq - ⟨basis, coords⟩ - -def vec (coords : List V) : NSpace coords.length V := coords.get - -------------------------------------------------------------------------------- def main : IO Unit :=