From 39bedbfd82f4ad5a07b9f91c5514f9a05d56d1e7 Mon Sep 17 00:00:00 2001 From: Michael Reed <18372368+chakravala@users.noreply.github.com> Date: Sun, 15 Dec 2019 17:07:06 -0500 Subject: [PATCH] fixed chain contraction product of tangent algebra --- .travis.yml | 2 +- README.md | 14 +- docs/make.jl | 4 +- docs/src/algebra.md | 251 +++++++------------- docs/src/design.md | 62 ++--- docs/src/index.md | 4 +- docs/src/library.md | 246 +------------------- docs/src/tutorials/algebra-of-space.md | 207 ++++++----------- docs/src/tutorials/install.md | 6 +- docs/src/tutorials/mixed-tensors.md | 43 ++-- docs/src/tutorials/quick-start.md | 42 ++-- src/algebra.jl | 305 +++++++++++++++++-------- 12 files changed, 438 insertions(+), 748 deletions(-) diff --git a/.travis.yml b/.travis.yml index 2996650..5f00399 100644 --- a/.travis.yml +++ b/.travis.yml @@ -20,6 +20,6 @@ jobs: julia: 1.3 os: linux script: - - julia --project=docs/ -e 'using Pkg; Pkg.add("Documenter"); Pkg.add("DirectSum"); Pkg.add("AbstractTensors"); Pkg.develop(PackageSpec(path=pwd()))' + - julia --project=docs/ -e 'using Pkg; Pkg.add("Documenter"); Pkg.add("DirectSum"); Pkg.add("AbstractTensors"); Pkg.add("Leibniz"); Pkg.develop(PackageSpec(path=pwd()))' - julia --project=docs/ docs/make.jl after_success: skip diff --git a/README.md b/README.md index d58794b..bb25ea2 100644 --- a/README.md +++ b/README.md @@ -74,7 +74,19 @@ More information and tutorials are available at https://grassmann.crucialflow.co ### Requirements -Availability of this package and its subpackages can be automatically handled with the Julia package manager; however, when the `master` branch is used it is possible that some of the dependencies also require a development branch before the release. This may include (but is not limited to) the following packages: +*Grassmann.jl* is a package for the [Julia language](https://julialang.org), which can be obtained from their website or the recommended method for your operating system (GNU/Linux/Mac/Windows). Go to [docs.julialang.org](https://docs.julialang.org) for documentation. + +Availability of this package and its subpackages can be automatically handled with the Julia package manager `using Pkg` and `Pkg.add("Grassmann")` or by entering: +```Julia +pkg> add Grassmann +``` +If you would like to keep up to date with the latest commits, instead use +```Julia +pkg> add Grassmann#master +``` +which is not recommended if you want to use a stable release. + +When the `master` branch is used it is possible that some of the dependencies also require a development branch before the release. This may include (but is not limited to) the following packages: This requires a merged version of `ComputedFieldTypes` at https://github.com/vtjnash/ComputedFieldTypes.jl diff --git a/docs/make.jl b/docs/make.jl index db3463d..e9d804d 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,11 +1,11 @@ # This file is part of Grassmann.jl. It is licensed under the GPL license # Grassmann Copyright (C) 2019 Michael Reed -using Documenter, DirectSum, AbstractTensors, Grassmann +using Documenter, DirectSum, AbstractTensors, Leibniz, Grassmann makedocs( # options - modules = [DirectSum,AbstractTensors,Grassmann], + modules = [DirectSum,AbstractTensors,Leibniz,Grassmann], doctest = false, format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"), sitename = "Grassmann.jl", diff --git a/docs/src/algebra.md b/docs/src/algebra.md index 88e2a0f..0d3cd0b 100644 --- a/docs/src/algebra.md +++ b/docs/src/algebra.md @@ -23,62 +23,33 @@ The formal sum of equal `grade` elements is an oriented `Chain` and with mixed ` Thus, various standard operations on the oriented multi-sets are possible including `∪,∩,⊕` and the index operation `⊖`, which is symmetric difference operation `⊻`. The elements of the `Algebra` can be generated in many ways using the `Basis` elements created by the `@basis` macro, -```Julia -julia> using Grassmann; @basis ℝ'⊕ℝ^3 # equivalent to basis"-+++" -(⟨-+++⟩, v, v₁, v₂, v₃, v₄, v₁₂, v₁₃, v₁₄, v₂₃, v₂₄, v₃₄, v₁₂₃, v₁₂₄, v₁₃₄, v₂₃₄, v₁₂₃₄) +```@repl ga +using Grassmann; @basis ℝ'⊕ℝ^3 # equivalent to basis"-+++" ``` As a result of this macro, all of the `Basis{V,G}` elements generated by that `VectorBundle` become available in the local workspace with the specified naming. The first argument provides signature specifications, the second argument is the variable name for the `VectorBundle`, and the third and fourth argument are prefixes of the `Basis` vector names (and covector basis names). By default, `V` is assigned the `VectorBundle` and `v` is the prefix for the `Basis` elements. -```Julia -julia> V # Minkowski spacetime -⟨-+++⟩ - -julia> typeof(V) # dispatch by vector space -VectorBundle{4,0,0x0000000000000001} - -julia> typeof(v13) # extensive type info -Basis{⟨-+++⟩,2,0x0000000000000005} - -julia> v13∧v2 # exterior tensor product --1v₁₂₃ - -julia> ans^2 # applies geometric product -1v - -julia> @btime 2v1+v3 # vector element - 37.794 ns (3 allocations: 80 bytes) -2v₁ + 0v₂ + 1v₃ + 0v₄ - -julia> @btime $ans⋅$ans # inner product - 15.266 ns (3 allocations: 48 bytes) --3v +```@repl ga +V # Minkowski spacetime +typeof(V) # dispatch by vector space +typeof(v13) # extensive type info +v13∧v2 # exterior tensor product +ans^2 # applies geometric product +2v1+v3 # vector element +ans⋅ans # inner product ``` It is entirely possible to assign multiple different bases with different signatures without any problems. In the following command, the `@basis` macro arguments are used to assign the vector space name to `S` instead of `V` and basis elements to `b` instead of `v`, so that their local names do not interfere: -```Julia -julia> @basis "++++" S b; - -julia> let k = (b1+b2)-b3 - for j ∈ 1:9 - k = k*(b234+b134) - println(k) - end end -0 + 1v₁₄ + 1v₂₄ + 2v₃₄ -0 - 2v₁ - 2v₂ + 2v₃ -0 - 2v₁₄ - 2v₂₄ - 4v₃₄ -0 + 4v₁ + 4v₂ - 4v₃ -0 + 4v₁₄ + 4v₂₄ + 8v₃₄ -0 - 8v₁ - 8v₂ + 8v₃ -0 - 8v₁₄ - 8v₂₄ - 16v₃₄ -0 + 16v₁ + 16v₂ - 16v₃ -0 + 16v₁₄ + 16v₂₄ + 32v₃₄ +```@repl ga +@basis "++++" S b; +let k = (b1+b2)-b3 + for j ∈ 1:9 + k = k*(b234+b134) + println(k) +end end ``` Alternatively, if you do not wish to assign these variables to your local workspace, the versatile `Grassmann.Algebra{N}` constructors can be used to contain them, which is exported to the user as the method `Λ(V)`, -```Julia -julia> G3 = Λ(3) # equivalent to Λ(V"+++"), Λ(ℝ^3), Λ.V3 -Grassmann.Algebra{⟨+++⟩,8}(v, v₁, v₂, v₃, v₁₂, v₁₃, v₂₃, v₁₂₃) - -julia> G3.v13 * G3.v12 -v₂₃ +```@repl ga +G3 = Λ(3) # equivalent to Λ(V"+++"), Λ(ℝ^3), Λ.V3 +G3.v13 * G3.v12 ``` The *geometric algebraic product* is the oriented symmetric difference operator `⊖` (weighted by the bilinear form `g`) and multi-set sum `⊕` applied to multilinear tensor products `⊗` in a single operation. Symmetry properties of the tensor algebra can be characterized in terms of the geometric product by two averaging operations, which are the symmetrization `⊙` and anti-symmetrization `⊠` operators. @@ -86,51 +57,28 @@ These products satisfy various `MultiVector` properties, including the associati It is possible to assign the **quaternion** generators `i,j,k` with -```Julia -julia> i,j,k = hyperplanes(ℝ^3) -3-element Array{Simplex{⟨+++⟩,2,B,Int64} where B,1}: - -1v₂₃ - 1v₁₃ - -1v₁₂ - -julia> @btime $i^2, $j^2, $k^2, $i*$j*$k - 0.027 ns (0 allocations: 0 bytes) -(-1v, -1v, -1v, -1v) - -julia> @btime -(j+k) * (j+k) - 97.373 ns (4 allocations: 176 bytes) -2.0v⃖ - -julia> @btime -(j+k) * i - 67.695 ns (3 allocations: 144 bytes) -0.0 - 1.0v₁₂ - 1.0v₁₃ +```@repl ga +i,j,k = hyperplanes(ℝ^3) +i^2, j^2, k^2, i*j*k +-(j+k) * (j+k) +-(j+k) * i ``` Alternatively, another representation of the quaternions is -```Julia -julia> basis"--" -(⟨--⟩, v, v₁, v₂, v₁₂) - -julia> v1^2, v2^2, v12^2, v1*v2*v12 -(-1v, -1v, -1v, -1v) +```@repl ga +basis"--" +v1^2, v2^2, v12^2, v1*v2*v12 ``` With the preliminary implementations of the `exp` and `log` functions one can compute -```Julia -julia> exp(0.5π/2*(i+j)/sqrt(2)) -0.7071067811865476 + 0.5v₁₃ - 0.5v₂₃ - -julia> ans == (sqrt(2)+j+i)/2 -true - -julia> log1p(i) -0.34657359027997264 - 0.7853981633974485v₂₃ - -julia> log(i) -0.0 - 1.5708963467978978v₂₃ +```@repl ga +exp(0.5π/2*(i+j)/sqrt(2)) +ans == (sqrt(2)+j+i)/2 +log1p(i) +log(i) ``` The parametric type formalism in `Grassmann` is highly expressive to enable the pre-allocation of geometric algebra computations for specific sparse-subalgebras, including the representation of rotational groups, Lie bivector algebras, and affine projective geometry. With [LightGraphs,jl](https://github.com/JuliaGraphs/LightGraphs.jl), [GraphPlot.jl](https://github.com/JuliaGraphs/GraphPlot.jl), [Cairo.jl](https://github.com/JuliaGraphics/Cairo.jl), [Compose.jl](https://github.com/GiovineItalia/Compose.jl) it is possible to convert `Grassmann` numbers into graphs. -```Julia +```julia using Grassmann, Compose # environment: LightGraphs, GraphPlot x = Grassmann.Algebra(ℝ^7).v123 Grassmann.graph(x+!x) @@ -139,7 +87,7 @@ draw(PDF("simplex.pdf",16cm,16cm),x+!x) ![paper/img/triangle-tetrahedron.png](https://raw.githubusercontent.com/chakravala/Grassmann.jl/master/paper/img/triangle-tetrahedron.png) Due to [GeometryTypes,jl](https://github.com/JuliaGeometry/GeometryTypes.jl) `Point` interoperability, plotting and visualizing with [Makie.jl](https://github.com/JuliaPlots/Makie.jl) is easily possible. For example, the `vectorfield` method creates an anonymous `Point` function that applies a versor outermorphism: -```Julia +```julia using Grassmann, Makie basis"2" # Euclidean streamplot(vectorfield(exp(π*v12/2)),-1.5..1.5,-1.5..1.5) @@ -154,7 +102,7 @@ streamplot(vectorfield(v1*exp((π/4)*v12/2)),-1.5..1.5,-1.5..1.5) ![paper/img/plane-3.png](https://raw.githubusercontent.com/chakravala/Grassmann.jl/master/paper/img/plane-3.png) ![paper/img/plane-4.png](https://raw.githubusercontent.com/chakravala/Grassmann.jl/master/paper/img/plane-4.png) ![paper/img/plane-3.png](https://raw.githubusercontent.com/chakravala/Grassmann.jl/master/paper/img/plane-5.png) ![paper/img/plane-4.png](https://raw.githubusercontent.com/chakravala/Grassmann.jl/master/paper/img/plane-6.png) -```Julia +```julia using Grassmann, Makie @basis S"∞+++" f(t) = (↓(exp(π*t*((3/7)*v12+v∞3))>>>↑(v1+v2+v3))) @@ -165,26 +113,26 @@ lines(points(f,V(3,4,5))) ``` ![paper/img/torus.png](https://raw.githubusercontent.com/chakravala/Grassmann.jl/master/paper/img/torus.png) ![paper/img/helix.png](https://raw.githubusercontent.com/chakravala/Grassmann.jl/master/paper/img/helix.png) -```Julia +```julia using Grassmann, Makie; @basis S"∞+++" streamplot(vectorfield(exp((π/4)*(v12+v∞3)),V(2,3,4)),-1.5..1.5,-1.5..1.5,-1.5..1.5,gridsize=(10,10)) ``` ![paper/img/orb.png](https://raw.githubusercontent.com/chakravala/Grassmann.jl/master/paper/img/orb.png) -```Julia +```julia using Grassmann, Makie; @basis S"∞+++" streamplot(vectorfield(exp((π/4)*(v12+v∞3)),V(2,3,4),V(1,2,3)),-1.5..1.5,-1.5..1.5,-1.5..1.5,gridsize=(10,10)) ``` ![paper/img/wave.png](https://raw.githubusercontent.com/chakravala/Grassmann.jl/master/paper/img/wave.png) -```Julia +```julia using Grassmann, Makie; @basis S"∞+++" f(t) = ↓(exp(t*v∞*(sin(3t)*3v1+cos(2t)*7v2-sin(5t)*4v3)/2)>>>↑(v1+v2-v3)) lines(points(f,V(2,3,4))) ``` ![paper/img/orb.png](https://raw.githubusercontent.com/chakravala/Grassmann.jl/master/paper/img/orbit-2.png) -```Julia +```julia using Grassmann, Makie; @basis S"∞+++" f(t) = ↓(exp(t*(v12+0.07v∞*(sin(3t)*3v1+cos(2t)*7v2-sin(5t)*4v3)/2))>>>↑(v1+v2-v3)) lines(points(f,V(2,3,4))) @@ -194,7 +142,7 @@ lines(points(f,V(2,3,4))) ## Approaching ∞ dimensions with `SparseAlgebra` and `ExtendedAlgebra` In order to work with a `TensorAlgebra{V}`, it is necessary for some computations to be cached. This is usually done automatically when accessed. -```Julia +```julia julia> Λ(7) + Λ(7)' Grassmann.SparseAlgebra{⟨+++++++-------⟩*,16384}(v, ..., v₁₂₃₄₅₆₇w¹²³⁴⁵⁶⁷) ``` @@ -204,21 +152,19 @@ Staging of precompilation and caching is designed so that a user can smoothly tr The parametric type formalism in `Grassmann` is highly expressive and enables pre-allocation of geometric algebra computations involving specific sparse subalgebras, including the representation of rotational groups. It is possible to reach `Simplex` elements with up to `N=62` vertices from a `TensorAlgebra` having higher maximum dimensions than supported by Julia natively. -```Julia -julia> Λ(62) -Grassmann.ExtendedAlgebra{⟨++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++⟩,4611686018427387904}(v, ..., v₁₂₃₄₅₆₇₈₉₀abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ) - -julia> Λ(62).v32a87Ng --1v₂₃₇₈agN +```@repl ga +Λ(62) +Λ(62).v32a87Ng ``` The 62 indices require full alpha-numeric labeling with lower-case and capital letters. This now allows you to reach up to `4,611,686,018,427,387,904` dimensions with Julia `using Grassmann`. Then the volume element is -```Julia -v₁₂₃₄₅₆₇₈₉₀abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ +```@example +using DirectSum # hide +DirectSum.printindices(stdout,DirectSum.indices(UInt(2^62-1))) # hide ``` Full `MultiVector` allocations are only possible for `N≤22`, but sparse operations are also available at higher dimensions. While `Grassmann.Algebra{V}` is a container for the `TensorAlgebra` generators of `V`, the `Grassmann.Algebra` is only cached for `N≤8`. For the range of dimensions `8 Λ(22) Grassmann.SparseAlgebra{⟨++++++++++++++++++++++⟩,4194304}(v, ..., v₁₂₃₄₅₆₇₈₉₀abcdefghijkl) ``` @@ -226,12 +172,9 @@ This is the largest `SparseAlgebra` that can be generated with Julia, due to arr To reach higher dimensions with `N>22`, the `Grassmann.ExtendedAlgebra` type is used. It is suficient to work with a 64-bit representation (which is the default). And it turns out that with 62 standard keyboard characters, this fits nicely. -```Julia -julia> V = ℝ^22 -⟨++++++++++++++++++++++⟩ - -julia> Λ(V+V') -Grassmann.ExtendedAlgebra{⟨++++++++++++++++++++++----------------------⟩*,17592186044416}(v, ..., v₁₂₃₄₅₆₇₈₉₀abcdefghijklw¹²³⁴⁵⁶⁷⁸⁹⁰ABCDEFGHIJKL) +```@repl ga +V = ℝ^22 +Λ(V+V') ``` At 22 dimensions and lower there is better caching, with further extra caching for 8 dimensions or less. Thus, the largest Hilbert space that is fully reachable has 4,194,304 dimensions, but we can still reach out to 4,611,686,018,427,387,904 dimensions with the `ExtendedAlgebra` built in. @@ -243,38 +186,19 @@ The sparse representations are a work in progress to be improved with time. ## Null-basis of the projective split In the following example, the null-basis from the projective split is used: -```Julia -julia> using Grassmann; @basis S"∞∅++" -(⟨∞∅++⟩, v, v∞, v∅, v₁, v₂, v∞∅, v∞₁, v∞₂, v∅₁, v∅₂, v₁₂, v∞∅₁, v∞∅₂, v∞₁₂, v∅₁₂, v∞∅₁₂) - -julia> v∞^2, v∅^2, v1^2, v2^2 -(0v, 0v, v, v) - -julia> v∞ ⋅ v∅ --1v - -julia> v∞∅^2 -v - -julia> v∞∅ * v∞, v∞∅ * v∅ -(-1v∞, v∅) - -julia> v∞ * v∅, v∅ * v∞ -(-1 + 1v∞∅, -1 - 1v∞∅) +```@repl ga +using Grassmann; @basis S"∞∅++" +v∞^2, v∅^2, v1^2, v2^2 +v∞ ⋅ v∅ +v∞∅^2 +v∞∅ * v∞, v∞∅ * v∅ +v∞ * v∅, v∅ * v∞ ``` -## Differential forms and Leibniz tangent algebra - -Multiplication with an `ϵᵢ` element is used help signify tensor fields so that differential operators are automatically applied in the `Basis` algebra as ∂ⱼ⊖(ω⊗ϵᵢ) = ∂ⱼ(ωϵᵢ) ≠ (∂ⱼ⊗ω)⊖ϵᵢ. -```Julia -julia> using Reduce, Grassmann; @mixedbasis tangent(ℝ^2,3,2); - -julia> (∂1+∂12) * (:(x1^2*x2^2)*ϵ1 + :(sin(x1))*ϵ2) -0.0 + (2 * x1 * x2 ^ 2)∂₁ϵ¹ + (cos(x1))∂₁ϵ² + (4 * x1 * x2)∂₁₂ϵ¹ -``` +## Differential forms and tangent algebra The product rule is encoded into `Grassmann` algebra when a `tangent` bundle is used, demonstrated here symbolically with `Reduce` by using the dual number definition: -```Julia +```julia julia> using Grassmann, Reduce Reduce (Free CSL version, revision 4590), 11-May-18 ... @@ -288,39 +212,30 @@ julia> a * b x * y + (dy * x + dx * y)v₁ϵ₁ ``` Higher order and multivariable Taylor numbers are also supported. -```Julia -julia> @basis tangent(ℝ,2,2) # 1D Grade, 2nd Order, 2 Variables -(⟨+₁₂⟩, v, v₁, ∂₁, ∂₂, ∂₁v₁, ∂₂v₁, ∂₁₂, ∂₁₂v₁) - -julia> ∂1 * ∂1v1 -∂₁∂₁v₁ - -julia> ∂1 * ∂2 -∂₁₂ - -julia> v1*∂12 -∂₁₂v₁ - -julia> ∂12*∂2 # 3rd order is zero -0v - -julia> @mixedbasis tangent(ℝ^2,2,2); # 2D Grade, 2nd Order, 2 Variables - -julia> ∇ = ∂1v1 + ∂2v2 # vector field -0v₁₂ + 1∂₁v₁ + 0∂₂v₁ + 0∂₁v₂ + 1∂₂v₂ + 0∂₁₂ - -julia> ∇ ⋅ ∇ # Laplacian -0.0v₁ + 0.0v₂ + 1∂₁∂₁ + 1∂₂∂₂ +```@repl ga +@basis tangent(ℝ,2,2) # 1D Grade, 2nd Order, 2 Variables +∂1 * ∂1v1 +∂1 * ∂2 +v1*∂12 +∂12*∂2 # 3rd order is zero +@mixedbasis tangent(ℝ^2,2,2); # 2D Grade, 2nd Order, 2 Variables +V(∇) # vector field +V(∇) ⋅ V(∇) # Laplacian +ans*∂1 # 3rd order is zero +``` +Multiplication with an `ϵᵢ` element is used help signify tensor fields so that differential operators are automatically applied in the `Basis` algebra as ∂ⱼ⊖(ω⊗ϵᵢ) = ∂ⱼ(ωϵᵢ) ≠ (∂ⱼ⊗ω)⊖ϵᵢ. +```julia +julia> using Reduce, Grassmann; @mixedbasis tangent(ℝ^2,3,2); -julia> ans*∂1 # 3rd order is zero -0.0v⃖ +julia> (∂1+∂12) * (:(x1^2*x2^2)*ϵ1 + :(sin(x1))*ϵ2) +0.0 + (2 * x1 * x2 ^ 2)∂₁ϵ¹ + (cos(x1))∂₁ϵ² + (4 * x1 * x2)∂₁₂ϵ¹ ``` Although fully generalized, the implementation in this release is still experimental. ## Symbolic coefficients by declaring an alternative scalar algebra Due to the abstract generality of the code generation of the `Grassmann` product algebra, it is easily possible to extend the entire set of operations to other kinds of scalar coefficient types. -```Julia +```julia julia> using GaloisFields, Grassmann julia> const F = GaloisField(7) @@ -329,22 +244,20 @@ julia> const F = GaloisField(7) julia> basis"2" (⟨++⟩, v, v₁, v₂, v₁₂) -julia> @btime F(3)*v1 - 21.076 ns (2 allocations: 32 bytes) +julia> F(3)*v1 3v₁ -julia> @btime inv($ans) - 14.965 ns (0 allocations: 0 bytes) +julia> inv(ans) 5v₁ ``` By default, the coefficients are required to be `<:Number`. However, if this does not suit your needs, alternative scalar product algebras can be specified with -```Julia +```julia Grassmann.generate_algebra(:AbstractAlgebra,:SetElem) ``` where `:SetElem` is the desired scalar field and `:AbstractAlgebra` is the scope which contains the scalar field. With the usage of `Requires`, symbolic scalar computation with [Reduce.jl](https://github.com/chakravala/Reduce.jl) and other packages is automatically enabled, -```Julia +```julia julia> using Reduce, Grassmann Reduce (Free CSL version, revision 4590), 11-May-18 ... @@ -362,7 +275,7 @@ a * c + b * d + (a * d - b * c)v₁₂ ``` If these compatibility steps are followed, then `Grassmann` will automatically declare the product algebra to use the `Reduce.Algebra` symbolic field operation scope. -```Julia +```julia julia> using Reduce,Grassmann; basis"4" Reduce (Free CSL version, revision 4590), 11-May-18 ... (⟨++++⟩, v, v₁, v₂, v₃, v₄, v₁₂, v₁₃, v₁₄, v₂₃, v₂₄, v₃₄, v₁₂₃, v₁₂₄, v₁₃₄, v₂₃₄, v₁₂₃₄) diff --git a/docs/src/design.md b/docs/src/design.md index d9c7bd0..f3f2671 100644 --- a/docs/src/design.md +++ b/docs/src/design.md @@ -20,7 +20,7 @@ Thus, computations involving fully general rotational algebras and Lie bivector Conformal geometric algebra is possible with the Minkowski plane ``v_{\infty\emptyset}``, based on the null-basis. In general, multivalued quantum logic is enabled by the ``\wedge,\vee,\star`` Grassmann lattice. Mixed-symmetry algebra with *Leibniz.jl* and *Grassmann.jl*, having the geometric algebraic product chain rule, yields automatic differentiation and Hodge-DeRahm co/homology as unveiled by Grassmann. -Most importantly, the Dirac-Clifford product yields generalized Hodge-Laplacian and the Betti numbers with Euler characteristic `χ`. +Most importantly, the Dirac-Clifford product yields generalized Hodge-Laplacian and the Betti numbers with Euler characteristic ``\chi``. Due to the abstract generality of the product algebra code generation, it is possible to extend the `Grassmann` library to include additional high performance products with few extra definitions. Operations on ultra-sparse representations for very high dimensional algebras will be gaining further performance enhancements in future updates, along with hybrid optimizations for low-dimensional algebra code generation. @@ -44,9 +44,11 @@ and `μ` is an integer specifying the order of the tangent bundle (i.e. multipli The metric signature of the `Basis{V,1}` elements of a vector space `V` can be specified with the `V"..."` constructor by using `+` and `-` to specify whether the `Basis{V,1}` element of the corresponding index squares to `+1` or `-1`. For example, `S"+++"` constructs a positive definite 3-dimensional `VectorBundle`. -```Julia -julia> ℝ^3 == V"+++" == vectorspace(3) -true +```@setup ds +using DirectSum +``` +```@repl ds +ℝ^3 == V"+++" == vectorspace(3) ``` It is also possible to specify an arbitrary `DiagonalForm` having numerical values for the basis with degeneracy `D"1,1,1,0"`, although the `Signature` format has a more compact representation. Further development will result in more metric types. @@ -55,49 +57,29 @@ Declaring an additional plane at infinity is done by specifying it in the string Additionally, the *null-basis* based on the projective split for confromal geometric algebra would be specified with `∞∅` initially (i.e. 5D CGA `S"∞∅+++"`). These two declared basis elements are interpreted in the type system. The `tangent` map takes `V` to its tangent space and can be applied repeatedly for higher orders, such that `tangent(V,μ,ν)` can be used to specify `μ` and `ν`. -```Julia -julia> V = tangent(ℝ^3) -⟨+++₁⟩ - -julia> V' -⟨---¹⟩' - -julia> V+V' -⟨+++---₁¹⟩* +```@repl ds +V = tangent(ℝ^3) +V' +V+V' ``` The direct sum operator `⊕` can be used to join spaces (alternatively `+`), and the dual space functor `'` is an involution which toggles a dual vector space with inverted signature. -```Julia -julia> V = ℝ'⊕ℝ^3 -⟨-+++⟩ - -julia> V' -⟨+---⟩' - -julia> W = V⊕V' -⟨-++++---⟩* +```@repl ds +V = ℝ'⊕ℝ^3 +V' +W = V⊕V' ``` The direct sum of a `VectorBundle` and its dual `V⊕V'` represents the full mother space `V*`. -```Julia -julia> collect(V) # all vector basis elements -Grassmann.Algebra{⟨-+++⟩,16}(v, v₁, v₂, v₃, v₄, v₁₂, v₁₃, v₁₄, v₂₃, v₂₄, v₃₄, v₁₂₃, v₁₂₄, v₁₃₄, ...) - -julia> collect(V') # all covector basis elements -Grassmann.Algebra{⟨+---⟩',16}(w, w¹, w², w³, w⁴, w¹², w¹³, w¹⁴, w²³, w²⁴, w³⁴, w¹²³, w¹²⁴, w¹³⁴, ...) - -julia> collect(W) # all mixed basis elements -Grassmann.Algebra{⟨-++++---⟩*,256}(v, v₁, v₂, v₃, v₄, w¹, w², w³, w⁴, v₁₂, v₁₃, v₁₄, v₁w¹, v₁w², ... +```@repl ds +collect(V) # all vector basis elements +collect(V') # all covector basis elements +collect(W) # all mixed basis elements ``` In addition to the direct-sum operation, several other operations are supported, such as `∪,∩,⊆,⊇` for set operations. Due to the design of the `VectorBundle` dispatch, these operations enable code optimizations at compile-time provided by the bit parameters. -```Julia -julia> ℝ+ℝ' ⊇ vectorspace(1) -true - -julia> ℝ ∩ ℝ' == vectorspace(0) -true - -julia> ℝ ∪ ℝ' == ℝ+ℝ' -true +```@repl ds +ℝ+ℝ' ⊇ vectorspace(1) +ℝ ∩ ℝ' == vectorspace(0) +ℝ ∪ ℝ' == ℝ+ℝ' ``` **Remark**. Although some of the operations like `∪` and `⊕` are similar and sometimes result in the same values, the `union` and `sum` are entirely different operations in general. diff --git a/docs/src/index.md b/docs/src/index.md index 093e37e..8c2786f 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -9,13 +9,13 @@ [![BiVector](https://img.shields.io/badge/bivector.net-Discourse-blueviolet)](https://bivector.net) The [Grassmann.jl](https://github.com/chakravala/Grassmann.jl) package provides tools for doing computations based on multi-linear algebra, differential geometry, and spin groups using the extended tensor algebra known as Leibniz-Grassmann-Clifford-Hestenes geometric algebra. -Combinatorial products include `∧, ∨, ⋅, *, ⋆, ', ~, d, ∂` (which are the exterior, regressive, inner, and geometric products; along with the Hodge star, adjoint, reversal, differential and boundary operators). +Combinatorial products include ``\wedge, \vee, \cdot, *, \star, ', \sim, d, \partial`` (which are the exterior, regressive, inner, and geometric products; along with the Hodge star, adjoint, reversal, differential and boundary operators). Kernelized operations are built up from composite sparse tensor products and Hodge duality, with high dimensional support for up to 62 indices using staged caching and precompilation. Code generation enables concise yet highly extensible definitions. [DirectSum.jl](https://github.com/chakravala/DirectSum.jl) multivector parametric type polymorphism is based on tangent bundle vector spaces and conformal projective geometry to make the dispatch highly extensible for many applications. Additionally, universal interoperability between different sub-algebras is enabled by [AbstractTensors.jl](https://github.com/chakravala/AbstractTensors.jl), on which the type system is built. ```@contents -Pages = ["index.md","design.md","algebra.md","library.md","references.md"] +Pages = ["design.md","algebra.md","library.md","references.md"] ``` This `Grassmann` package for the Julia language was created by [github.com/chakravala](https://github.com/chakravala) for mathematics and computer algebra research with differential geometric algebras. diff --git a/docs/src/library.md b/docs/src/library.md index b92e94f..f9b51cc 100644 --- a/docs/src/library.md +++ b/docs/src/library.md @@ -1,251 +1,9 @@ # Grassmann.jl Library -```@contents -Pages = ["index.md","design.md","algebra.md","library.md","references.md"] -``` - -## Index - ```@index Pages = ["library.md"] ``` -## DirectSum Interface - -```@docs -DirectSum.Manifold -``` - -```@docs -DirectSum.VectorBundle -``` - -```@docs -DirectSum.Signature -``` - -```@docs -DirectSum.DiagonalForm -``` - -```@docs -DirectSum.SubManifold -``` - -```@docs -DirectSum.vectorspace -``` - -```@docs -DirectSum.diagonalform +```@autodocs +Modules = [DirectSum, Grassmann, AbstractTensors, Leibniz] ``` - -```@docs -DirectSum.var"@V_str" -``` - -```@docs -DirectSum.var"@S_str" -``` - -```@docs -DirectSum.var"@D_str" -``` - -```@docs -DirectSum.value -``` - -```@docs -DirectSum.grade -``` - -```@docs -DirectSum.tangent -``` - -```@docs -DirectSum.indices -``` - -```@docs -DirectSum.printindices -``` - -```@docs -DirectSum.printlabel -``` - -```@docs -DirectSum.indexstring -``` - -```@docs -DirectSum.indexsymbol -``` - -```@docs -DirectSum.:⊕ -``` - -## AbstractTensors Interface - -```@docs -AbstractTensors.TensorAlgebra -``` - -```@docs -AbstractTensors.interop -``` - -```@docs -AbstractTensors.interform -``` - -## Grassmann Interface - -```@docs -Grassmann.TensorGraded -``` - -```@docs -Grassmann.TensorTerm -``` - -```@docs -Grassmann.TensorMixed -``` - -```@docs -Grassmann.Basis -``` - -```@docs -Grassmann.Simplex -``` - -```@docs -Grassmann.Chain -``` - -```@docs -Grassmann.MultiVector -``` - -```@docs -Grasmann.SparseChain -``` - -```@docs -Grassmann.MultiGrade -``` - -```@docs -Grassmann.valuetype -``` - -```@docs -Grassmann.basis -``` - -```@docs -Grassmann.order -``` - -```@docs -Grassmann.χ -``` - -```@docs -Grassmann.gdims -``` - -```@docs -Grassmann.betti -``` - -```@docs -Grassmann.scalar -``` - -```@docs -Grassmann.vector -``` - -```@docs -Grassmann.volume -``` - -```@docs -Grassmann.isscalar -``` - -```@docs -Grassmann.angular -``` - -```@docs -Grassmann.radial -``` - -```@docs -Grassmann.complementrighthodge -``` - -```@docs -Grassmann.complementlefthodge -``` - -```@docs -Grassmann.complementright -``` - -```@docs -Grassmann.complementleft -``` - -```@docs -Grassmann.reverse -``` - -```@docs -Grassmann.involute -``` - -```@docs -Grassmann.:∗ -``` - -```@docs -Grassmann.:⊛ -``` - -```@docs -Grassmann.:∧ -``` - -```@docs -Grassmann.:∨ -``` - -```@docs -Grassmann.contraction -``` - -```@docs -Grassmann.:⊙ -``` - -```@docs -Grassmann.:⊠ -``` - -```@docs -Grassmann.:⊘ -``` - -```@docs -Grassmann.generate_algebra -``` - -etc, ... diff --git a/docs/src/tutorials/algebra-of-space.md b/docs/src/tutorials/algebra-of-space.md index c014755..6268fa4 100644 --- a/docs/src/tutorials/algebra-of-space.md +++ b/docs/src/tutorials/algebra-of-space.md @@ -4,138 +4,92 @@ This notebook is an adaptation from the [clifford](https://clifford.readthedocs. Import `Grassmann` and instantiate a three dimensional geometric algebra -```Julia -julia> using Grassmann - -julia> basis"3" -(⟨+++⟩, v, v₁, v₂, v₃, v₁₂, v₁₃, v₂₃, v₁₂₃) +```@repl ga +using Grassmann +basis"3" ``` Given a three dimensional GA with the orthonormal basis ``v_i\cdot v_j = \delta_{ij}``, the basis consists of scalars, three vectors, three bivectors, and a trivector. ```math -\{\underbrace{v}{\mbox{scalar}},\qquad\underbrace{v_1,v_2,v_3}{\mbox{vectors}},\qquad\underbrace{v_{12},v_{23},v_{13}}{\mbox{bivectors}},\qquad\underbrace{v_{123}}{\mbox{trivector}}\} +\{\underbrace{v}_{\text{scalar}},\qquad\underbrace{v_1,v_2,v_3}_{\text{vectors}},\qquad\underbrace{v_{12},v_{23},v_{13}}_{\text{bivectors}},\qquad\underbrace{v_{123}}_{\text{trivector}}\} ``` The `@basis` macro declares the algebra and assigns the `Basis` elements to local variables. The `Grassmann.Algebra` can also be assigned to `G3` as -```Julia -julia> G3 = Λ(3) -Grassmann.Algebra{⟨+++⟩,8}(v, v₁, v₂, v₃, v₁₂, v₁₃, v₂₃, v₁₂₃) +```@repl ga +G3 = Λ(3) ``` You may wish to explicitly assign the blades to variables like so, -```Julia +```julia e1 = G3.v1 e2 = G3.v2 # etc ... ``` Or, if you're lazy you can use the macro with different local names -```Julia -julia> @basis ℝ^3 E e -(⟨+++⟩, v, v₁, v₂, v₃, v₁₂, v₁₃, v₂₃, v₁₂₃) - -julia> e3, e123 -(v₃, v₁₂₃) +```@repl ga +@basis ℝ^3 E e +e3, e123 ``` ## Basics The basic products are available -```Julia -julia> v1*v2 # geometric product -v₁₂ - -julia> v1|v2 # inner product -0v - -julia> v1 ∧ v2 # exterior product -v₁₂ - -julia> v1 ∧ v2 ∧ v3 # even more exterior products -v₁₂₃ +```@repl ga +v1 * v2 # geometric product +v1 | v2 # inner product +v1 ∧ v2 # exterior product +v1 ∧ v2 ∧ v3 # even more exterior products ``` Multivectors can be defined in terms of the basis blades. For example, you can construct a rotor as a sum of a scalar and a bivector, like so -```Julia -julia> θ = π/4 -0.7853981633974483 - -julia> R = cos(θ) - sin(θ)*v23 -0.7071067811865476 - 0.7071067811865475v₂₃ - -julia> R = ~exp(θ*v23) -0.7071067811865476 - 0.7071067811865475v₂₃ +```@repl ga +θ = π/4 +R = cos(θ) + sin(θ)*v23 +R = exp(θ*v23) ``` You can also mix grades without any reason -```Julia -julia> A = 1 + 2v1 + 3v12 + 4v123 -1 + 2v₁ + 3v₁₂ + 4v₁₂₃ +```@repl ga +A = 1 + 2v1 + 3v12 + 4v123 ``` The reversion operator is accomplished with the tilde `~` in front of the `MultiVector` on which it acts -```Julia -julia> ~A -1 + 2v₁ - 3v₁₂ - 4v₁₂₃ +```@repl ga +~A ``` Taking a projection into a specific `grade` of a `MultiVector` is usually written ``\langle A\rangle_n`` and can be done using the soft brackets, like so -```Julia -julia> A(0) -1v - -julia> A(1) -2v₁ + 0v₂ + 0v₃ - -julia> A(2) -3v₁₂ + 0v₁₃ + 0v₂₃ +```@repl ga +A(0) +A(1) +A(2) ``` Using the reversion and grade projection operators, we can define the magnitude of `A` as ``|A|^2 = \langle\tilde A A\rangle`` -```Julia -julia> ~A*A -30 + 4v₁ + 12v₂ + 24v₃ - -julia> scalar(ans) -30v +```@repl ga +~A*A +scalar(ans) ``` This is done in the `abs` and `abs2` operators -```Julia -julia> abs2(A) -30 + 4v₁ + 12v₂ + 24v₃ - -julia> scalar(ans) -30v +```@repl ga +abs2(A) +scalar(ans) ``` The dual of a multivector `A` can be defined as ``\tilde AI``, where `I` is the pseudoscalar for the geometric algebra. In `G3`, the dual of a vector is a bivector: -```Julia -julia> a = 1v1 + 2v2 + 3v3 -1v₁ + 2v₂ + 3v₃ - -julia> ⋆a -3v₁₂ - 2v₁₃ + 1v₂₃ +```@repl ga +a = 1v1 + 2v2 + 3v3 +⋆a ``` ## Reflections Reflecting a vector ``c`` about a normalized vector ``n`` is pretty simple, ``c\mapsto -ncn`` -```Julia -julia> c = v1+v2+v3 # a vector -1v₁ + 1v₂ + 1v₃ - -julia> n = v1 # the reflector -v₁ - -julia> -n*c*n # reflect a in hyperplane normal to n -0.0 - 1.0v₁ + 1.0v₂ + 1.0v₃ +```@repl ga +c = v1+v2+v3 # a vector +n = v1 # the reflector +-n*c*n # reflect a in hyperplane normal to n ``` Because we have the `inv` available, we can equally well reflect in un-normalized vectors using ``a\mapsto n^{-1}an`` -```Julia -julia> a = v1+v2+v3 # the vector -1v₁ + 1v₂ + 1v₃ - -julia> n = 3v1 # the reflector -3v₁ - -julia> inv(n)*a*n -0.0 + 1.0v₁ - 1.0v₂ - 1.0v₃ - -julia> n\a*n -0.0 + 1.0v₁ - 1.0v₂ - 1.0v₃ +```@repl ga +a = v1+v2+v3 # the vector +n = 3v1 # the reflector +inv(n)*a*n +n\a*n ``` Reflections can also be made with respect to the hyperplane normal to the vector, in which case the formula is negated. @@ -143,72 +97,53 @@ Reflections can also be made with respect to the hyperplane normal to the vector A vector can be rotated using the formula ``a\mapsto \tilde R aR``, where `R` is a rotor. A rotor can be defined by multiple reflections, ``R = mn`` or by a plane and an angle ``R = e^{\theta B/2}``. For example, -```Julia -julia> R = exp(π/4*v12) -0.7071067811865476 + 0.7071067811865475v₁₂ - -julia> ~R*v1*R -0.0 + 2.220446049250313e-16v₁ + 1.0v₂ +```@repl ga +R = exp(π/4*v12) +~R*v1*R ``` Maybe we want to define a function which can return rotor of some angle ``\theta`` in the ``v_{12}``-plane, ``R_{12} = e^{\theta v_{12}/2}`` -```Julia +```@example ga R12(θ) = exp(θ/2*v12) +nothing # hide ``` And use it like this -```Julia -julia> R = R12(π/2) -0.7071067811865476 + 0.7071067811865475v₁₂ - -julia> a = v1+v2+v3 -1v₁ + 1v₂ + 1v₃ - -julia> ~R*a*R -0.0 - 0.9999999999999997v₁ + 1.0v₂ + 1.0v₃ +```@repl ga +R = R12(π/2) +a = v1+v2+v3 +~R*a*R ``` You might as well make the angle argument a bivector, so that you can control the plane of rotation as well as the angle -```Julia +```@example ga R_B(B) = exp(B/2) +nothing # hide ``` Then you could do -```Julia -julia> Rxy = R_B(π/4*v12) -0.9238795325112867 + 0.3826834323650898v₁₂ - -julia> Ryz = R_B(π/5*v23) -0.9510565162951535 + 0.3090169943749474v₂₃ +```@repl ga +Rxy = R_B(π/4*v12) +Ryz = R_B(π/5*v23) ``` or -```Julia -julia> R_B(π/6*(v23+v12)) -0.9322404424570728 + 0.25585909935689327v₁₂ + 0.25585909935689327v₂₃ +```@repl ga +R_B(π/6*(v23+v12)) ``` Maybe you want to define a function which returns a *function* that enacts a specified rotation, ``f(B) = a\mapsto e^{B/2}\\ae^{B/2}``. This just saves you having to write out the sandwich product, which is nice if you are cascading a bunch of rotors, like so -```Julia +```@example ga R_factory(B) = (R = exp(B/2); a -> ~R*a*R) Rxy = R_factory(π/3*v12) Ryz = R_factory(π/3*v23) Rxz = R_factory(π/3*v13) +nothing # hide ``` Then you can do things like -```Julia -julia> R = R_factory(π/6*(v23+v12)) # this returns a function -#7 (generic function with 1 method) - -julia> R(a) # which acts on a vector -0.0 + 0.5229556000177233v₁ + 0.7381444851051178v₂ + 1.4770443999822769v₃ + 2.7755575615628914e-17v₁₂₃ - -julia> Rxy(Ryz(Rxz(a))) -0.0 + 0.40849364905389035v₁ - 0.6584936490538903v₂ + 1.5490381056766584v₃ +```@repl ga +R = R_factory(π/6*(v23+v12)) # this returns a function +R(a) # which acts on a vector +Rxy(Ryz(Rxz(a))) ``` To make cascading a sequence of rotations as concise as possible, we could define a function which takes a list of bivectors ``A,B,C,...``, and enacts the sequence of rotations which they represent on some vector ``x``. -```Julia -julia> R_seq(args...) = (R = prod(exp.(args./2)); a -> ~R*a*R) -R_seq (generic function with 1 method) - -julia> R = R_seq(π/2*v23, π/2*v12, v1) -#11 (generic function with 1 method) - -julia> R(v1) -2.220446049250313e-16 + 3.469446951953614e-16v₁ + 0.9999999999999996v₂ - 1.3877787807814457e-17v₃ - 5.551115123125783e-17v₂₃ + 2.7755575615628914e-17v₁₂₃ +```@repl ga +R_seq(args...) = (R = prod(exp.(args./2)); a -> ~R*a*R) +R = R_seq(π/2*v23, π/2*v12, v1) +R(v1) ``` diff --git a/docs/src/tutorials/install.md b/docs/src/tutorials/install.md index 9ad7d50..c0112b1 100644 --- a/docs/src/tutorials/install.md +++ b/docs/src/tutorials/install.md @@ -1,8 +1,8 @@ # Installation -*Grassmann.jl* is a package for the [Julia language](https://julialang.org), which should be installed from website or the recommended method for your operating system (GNU/Linux/Mac/Windows). +*Grassmann.jl* is a package for the [Julia language](https://julialang.org), which can be obtained from their website or the recommended method for your operating system (GNU/Linux/Mac/Windows). Go to [docs.julialang.org](https://docs.julialang.org) for documentation. -Availability of this package and its subpackages can be automatically handled with the Julia package manager with `using Pkg` and `Pkg.add("Grassmann")` or by entering: +Availability of this package and its subpackages can be automatically handled with the Julia package manager `using Pkg` and `Pkg.add("Grassmann")` or by entering: ```Julia pkg> add Grassmann ``` @@ -18,7 +18,7 @@ When the `master` branch is used it is possible that some of the dependencies al This requires a merged version of `ComputedFieldTypes` at [ComputedFieldTypes.jl](https://github.com/vtjnash/ComputedFieldTypes.jl). -Interoperability of `TensorAlgebra` with other packages is automatically enabled by [DirectSum.jl](https://github.com/chakravala/DirectSum.jl) and [AbstractTensors.jl](https://github.com/chakravala/AbstractTensors.jl). +Interoperability of `TensorAlgebra` with other packages is enabled by [DirectSum.jl](https://github.com/chakravala/DirectSum.jl) and [AbstractTensors.jl](https://github.com/chakravala/AbstractTensors.jl). The package is compatible via [Requires.jl](https://github.com/MikeInnes/Requires.jl) with [Reduce.jl](https://github.com/chakravala/Reduce.jl), diff --git a/docs/src/tutorials/mixed-tensors.md b/docs/src/tutorials/mixed-tensors.md index df0b018..5bc2493 100644 --- a/docs/src/tutorials/mixed-tensors.md +++ b/docs/src/tutorials/mixed-tensors.md @@ -6,45 +6,34 @@ The package `Grassmann` is working towards making the full extent of this number ## Constructing linear transformations Note that `Λ(3)` gives the vector basis, and `Λ(3)'` gives the covector basis: -```Julia -julia> Λ(3) -Grassmann.Algebra{⟨+++⟩,8}(v, v₁, v₂, v₃, v₁₂, v₁₃, v₂₃, v₁₂₃) - -julia> Λ(3)' -Grassmann.Algebra{⟨---⟩',8}(w, w¹, w², w³, w¹², w¹³, w²³, w¹²³) +```@setup ga +using Grassmann +``` +```@repl ga +Λ(3) +Λ(3)' ``` The following command yields a local 2D vector and covector basis, -```Julia -julia> mixedbasis"2" -(⟨++--⟩*, v, v₁, v₂, w¹, w², v₁₂, v₁w¹, v₁w², v₂w¹, v₂w², w¹², v₁₂w¹, v₁₂w², v₁w¹², v₂w¹², v₁₂w¹²) - -julia> w1+2w2 -1w¹ + 2w² - -julia> ans(v1+v2) -3v +```@repl ga +mixedbasis"2" +w1+2w2 +ans(v1+v2) ``` The sum `w1+2w2` is interpreted as a covector element of the dual vector space, which can be evaluated as a linear functional when a vector argument is input. Using these in the workspace, it is possible to use the Grassmann exterior `∧`-tensor product operation to construct elements `ℒ` of the (1,1)-bivector subspace of linear transformations from the `Grade{2}` algebra. -```Julia -julia> ℒ = (v1+2v2)∧(3w1+4w2) -0v₁₂ + 3v₁w¹ + 4v₁w² + 6v₂w¹ + 8v₂w² + 0w¹² +```@repl ga +ℒ = (v1+2v2)∧(3w1+4w2) ``` The element `ℒ` is a linear form which can take `Grade{1}` vectors as input, -```Julia -julia> ℒ(v1+v2) -7v₁ + 14v₂ + 0w¹ + 0w² - -julia> L = [1,2] * [3,4]'; L * [1,1] -2-element Array{Int64,1}: - 7 - 14 +```@repl ga +ℒ(v1+v2) +L = [1,2] * [3,4]'; L * [1,1] ``` which is a computation equivalent to a matrix computation. The `TensorAlgebra` evalution is still a work in progress, and the API and implementation may change as more features and algebraic operations and product structure are added. -## Importing the generators of the Leech lattice +## Importing the Leech lattice generator In the example below, we define a constant `Leech` which can be used to obtain linear combinations of the Leech lattice, ```Julia diff --git a/docs/src/tutorials/quick-start.md b/docs/src/tutorials/quick-start.md index 18a3641..573ccbd 100644 --- a/docs/src/tutorials/quick-start.md +++ b/docs/src/tutorials/quick-start.md @@ -2,41 +2,25 @@ Import the `Grassmann` package and instantiate a two-dimensional algebra (G2), -```Julia -julia> using Grassmann - -julia> @basis ℝ^2 -(⟨++⟩, v, v₁, v₂, v₁₂) - -julia> v1*v2 # geometric product -v₁₂ - -julia> v1|v2 # inner product -0v - -julia> v1∧v2 # exterior product -v₁₂ +```@repl ga +using Grassmann +@basis ℝ^2 +v1*v2 # geometric product +v1|v2 # inner product +v1∧v2 # exterior product ``` ## Reflection -```Julia -julia> a = v1+v2 -1v₁ + 1v₂ - -julia> n = v1 -v₁ - -julia> -n*a/n # reflect a in hyperplane normal to n -0.0 - 1.0v₁ + 1.0v₂ +```@example ga +a = v1+v2 +n = v1 +-n*a/n # reflect a in hyperplane normal to n ``` ## Rotation -```Julia -julia> R = exp(π/4*v12) -0.7071067811865476 + 0.7071067811865475v₁₂ - -julia> R*v1*~R -0.0 + 2.220446049250313e-16v₁ - 1.0v₂ +```@repl ga +R = exp(π/4*v12) +~R*v1*R ``` diff --git a/src/algebra.jl b/src/algebra.jl index d27fc91..dddd1ab 100644 --- a/src/algebra.jl +++ b/src/algebra.jl @@ -292,13 +292,15 @@ function generate_mutators(M,F,set_val,SUB,MUL) end end -@inline exteraddmulti!(V::W,out,α,β,γ) where W<:Manifold = (count_ones(α&β)==0) && joinaddmulti!(V,out,α,β,γ) +@inline exterbits(V,α,β) = diffvars(V)≠0 ? ((a,b)=symmetricmask(V,α,β);count_ones(a&b)==0) : count_ones(α&β)==0 -@inline outeraddblade!(V::W,out,α,β,γ) where W<:Manifold = (count_ones(α&β)==0) && joinaddblade!(V,out,α,β,γ) +@inline exteraddmulti!(V::W,out,α,β,γ) where W<:Manifold = exterbits(V,α,β) && joinaddmulti!(V,out,α,β,γ) -@inline exteraddmulti!_pre(V::W,out,α,β,γ) where W<:Manifold = (count_ones(α&β)==0) && joinaddmulti!_pre(V,out,α,β,γ) +@inline outeraddblade!(V::W,out,α,β,γ) where W<:Manifold = exterbits(V,α,β) && joinaddblade!(V,out,α,β,γ) -@inline outeraddblade!_pre(V::W,out,α,β,γ) where W<:Manifold = (count_ones(α&β)==0) && joinaddblade!_pre(V,out,α,β,γ) +@inline exteraddmulti!_pre(V::W,out,α,β,γ) where W<:Manifold = exterbits(V,α,β) && joinaddmulti!_pre(V,out,α,β,γ) + +@inline outeraddblade!_pre(V::W,out,α,β,γ) where W<:Manifold = exterbits(V,α,β) && joinaddblade!_pre(V,out,α,β,γ) # Hodge star ★ @@ -793,273 +795,388 @@ function generate_products(Field=Field,VEC=:mvec,MUL=:*,ADD=:+,SUB=:-,CONJ=:conj #∧(a::$Field,b::Chain{V,G,T}) where {V,G,T<:$Field} = Chain{V,G,T}(a.*b.v) #∧(a::Chain{V,G,T},b::$Field) where {V,G,T<:$Field} = Chain{V,G,T}(a.v.*b) @generated function contraction(a::Chain{V,G,T},b::Basis{V,L}) where {V,G,T<:$Field,L} - Gndims(V) && (return g_zero(V)) + G+L>ndims(V) && diffvars(V)==0 && (return g_zero(V)) if binomial(ndims(w),G)<(1<0,mixedmode(W)>0 ? dual(V,bits(b)) : bits(b) for i ∈ 1:binomial(ndims(w),G) X = @inbounds C ? dual(V,ib[i]) : ib[i] - @inbounds outeraddblade!_pre(V,out,X,y,derive_pre(V,X,y,:(a[$i]),true)) + if μ + exteraddmulti!_pre(V,out,X,y,derive_pre(V,X,y,:(a[$i]),true)) + else + outeraddblade!_pre(V,out,X,y,derive_pre(V,X,y,:(a[$i]),true)) + end + end + return if μ + :(MultiVector{$V}($(Expr(:call,:SVector,out...)))) + else + :(Chain{$V,G+L,promote_type(T,valuetype(b))}($(Expr(:call,:SVector,out...)))) end - return :(Chain{$V,G+L,promote_type(T,valuetype(b))}($(Expr(:call,:SVector,out...)))) else return quote V = $V $(insert_expr((:N,:t,:μ),$(QuoteNode(VEC)))...) ib = indexbasis(ndims(w),G) - out = zeros($$VEC(N,G+L,t)) + out = zeros(μ ? $$VEC(N,t) : $$VEC(N,G+L,t)) C,y = mixedmode(w)>0,mixedmode(W)>0 ? dual(V,bits(b)) : bits(b) for i ∈ 1:binomial(ndims(w),G) X = @inbounds C ? dual(V,ib[i]) : ib[i] - if @inbounds outeraddblade!(V,out,X,y,derive_mul(V,X,y,a[i],true))&μ - out,t = zeros(svec(N,G+L,Any)) .+ out,Any + if μ + if @inbounds exteraddmulti!(V,out,X,y,derive_mul(V,X,y,a[i],true)) + out,t = zeros(svec(N,Any)) .+ out,Any + @inbounds exteraddmulti!(V,out,X,y,derive_mul(V,X,y,a[i],true)) + end + else @inbounds outeraddblade!(V,out,X,y,derive_mul(V,X,y,a[i],true)) end end - return Chain{V,G+L,t}(out) + return μ ? MultiVector{V}(out) : Chain{V,G+L,t}(out) end end end @generated function ∧(a::Basis{w,G},b::Chain{W,L,T}) where {w,W,T<:$Field,G,L} V = w==W ? w : ((w==dual(W)) ? (mixedmode(w)≠0 ? W+w : w+W) : (return :(interop(∧,a,b)))) - G+L>ndims(V) && (return g_zero(V)) + G+L>ndims(V) && diffvars(V)==0 && (return g_zero(V)) if binomial(ndims(W),L)<(1<0,mixedmode(w)>0 ? dual(V,bits(a)) : bits(a) for i ∈ 1:binomial(ndims(W),L) X = @inbounds C ? dual(V,ib[i]) : ib[i] - @inbounds outeraddblade!_pre(V,out,x,X,derive_pre(V,x,X,:(b[$i]),false)) + if μ + exteraddmulti!_pre(V,out,x,X,derive_pre(V,x,X,:(b[$i]),false)) + else + outeraddblade!_pre(V,out,x,X,derive_pre(V,x,X,:(b[$i]),false)) + end + end + return if μ + :(MultiVector{$V}($(Expr(:call,:SVector,out...)))) + else + :(Chain{$V,G+L,promote_type(T,valuetype(a))}($(Expr(:call,:SVector,out...)))) end - return :(Chain{$V,G+L,promote_type(T,valuetype(a))}($(Expr(:call,:SVector,out...)))) else return quote V = $V $(insert_expr((:N,:t,:μ),$(QuoteNode(VEC)))...) ib = indexbasis(ndims(W),L) - out = zeros($$VEC(N,G+L,t)) + out = zeros(μ ? $$VEC(N,t) : $$VEC(N,G+L,t)) C,x = mixedmode(W)>0,mixedmode(w)>0 ? dual(V,bits(a)) : bits(a) for i ∈ 1:binomial(ndims(W),L) X = @inbounds C ? dual(V,ib[i]) : ib[i] - if @inbounds outeraddblade!(V,out,x,X,derive_mul(V,x,X,b[i],false))&μ - out,t = zeros(svec(N,G+L,Any)) .+ out,Any + if μ + if @inbounds exteraddmulti!(V,out,x,X,derive_mul(V,x,X,b[i],false)) + out,t = zeros(svec(N,Any)) .+ out,Any + @inbounds eteraddmulti!(V,out,x,X,derive_mul(V,x,X,b[i],false)) + end + else @inbounds outeraddblade!(V,out,x,X,derive_mul(V,x,X,b[i],false)) end end - return Chain{V,G+L,t}(out) + return μ ? MultiVector{V}(out) : Chain{V,G+L,t}(out) end end end @generated function contraction(a::Chain{V,G,T},b::Simplex{V,L,B,S}) where {V,G,T<:$Field,B,S<:$Field,L} - Gndims(V) && (return g_zero(V)) + G+L>ndims(V) && diffvars(V)==0 && (return g_zero(V)) if binomial(ndims(w),G)<(1<0,mixedmode(W)>0 ? dual(V,bits(B)) : bits(B) for i ∈ 1:binomial(ndims(w),G) X = @inbounds C ? dual(V,ib[i]) : ib[i] - @inbounds outeraddblade!_pre(V,out,X,y,derive_pre(V,X,y,:(a[$i]),:(b.v),$(QuoteNode(MUL)))) + if μ + exteraddmulti!_pre(V,out,X,y,derive_pre(V,X,y,:(a[$i]),:(b.v),$(QuoteNode(MUL)))) + else + outeraddblade!_pre(V,out,X,y,derive_pre(V,X,y,:(a[$i]),:(b.v),$(QuoteNode(MUL)))) + end + end + return if μ + :(multivector{$v}($(expr(:call,:svector,out...)))) + else + :(Chain{$V,G+L,promote_type(T,S)}($(Expr(:call,:SVector,out...)))) end - return :(Chain{$V,G+L,promote_type(T,S)}($(Expr(:call,:SVector,out...)))) else return quote $(insert_expr((:N,:t,:μ),$(QuoteNode(VEC)))...) ib = indexbasis(ndims(w),G) - out = zeros($$VEC(N,G+L,t)) + out = zeros(μ ? $$VEC(N,t) : $$VEC(N,G+L,t)) C,y = mixedmode(w)>0,mixedmode(W)>0 ? dual(V,bits(B)) : bits(B) for i ∈ 1:binomial(ndims(w),G) X = @inbounds C ? dual(V,ib[i]) : ib[i] - if @inbounds outeraddblade!(V,out,X,y,derive_mul(V,X,y,a[i],b.v,$$MUL))&μ - out,t = zeros(svec(N,G+L,Any)) .+ out,Any + if μ + if @inbounds exteraddmulti!(V,out,X,y,derive_mul(V,X,y,a[i],b.v,$$MUL)) + out,t = zeros(svec(N,Any)) .+ out,Any + @inbounds exteraddmulti!(V,out,X,y,derive_mul(V,X,y,a[i],b.v,$$MUL)) + end + else @inbounds outeraddblade!(V,out,X,y,derive_mul(V,X,y,a[i],b.v,$$MUL)) end end - return Chain{V,G+L,t}(out) + return μ ? MultiVector{V}(out) : Chain{V,G+L,t}(out) end end end @generated function ∧(a::Simplex{w,G,B,S},b::Chain{W,L,T}) where {T<:$Field,w,W,B,S<:$Field,G,L} V = w==W ? w : ((w==dual(W)) ? (mixedmode(w)≠0 ? W+w : w+W) : (return :(interop(∧,a,b)))) - G+L>ndims(V) && (return g_zero(V)) + G+L>ndims(V) && diffvars(V)==0 && (return g_zero(V)) if binomial(ndims(W),L)<(1<0,mixedmode(w)>0 ? dual(V,bits(B)) : bits(B) for i ∈ 1:binomial(ndims(W),L) X = @inbounds C ? dual(V,ib[i]) : ib[i] - @inbounds outeraddblade!_pre(V,out,x,X,derive_pre(V,x,X,:(a.v),:(b[$i]),$(QuoteNode(MUL)))) + if μ + exteraddmulti!_pre(V,out,x,X,derive_pre(V,x,X,:(a.v),:(b[$i]),$(QuoteNode(MUL)))) + else + outeraddblade!_pre(V,out,x,X,derive_pre(V,x,X,:(a.v),:(b[$i]),$(QuoteNode(MUL)))) + end + end + return if μ + :(MultiVector{$V}($(Expr(:call,:SVector,out...)))) + else + :(Chain{$V,G+L,promote_type(T,S)}($(Expr(:call,:SVector,out...)))) end - return :(Chain{$V,G+L,promote_type(T,S)}($(Expr(:call,:SVector,out...)))) else return quote $(insert_expr((:N,:t,:μ),$(QuoteNode(VEC)))...) ib = indexbasis(ndims(W),L) - out = zeros($$VEC(N,G+L,t)) + out = zeros(μ ? $$VEC(N,t) : $$VEC(N,G+L,t)) C,x = mixedmode(W)>0,mixedmode(w)>0 ? dual(V,bits(B)) : bits(B) for i ∈ 1:binomial(ndims(W),L) X = @inbounds C ? dual(V,ib[i]) : ib[i] - if @inbounds outeraddblade!(V,out,x,X,derive_mul(V,x,X,a.v,b[i],$$MUL))&μ - out,t = zeros(svec(N,G+L,Any)) .+ out,Any + if μ + if @inbounds exteraddmulti!(V,out,x,X,derive_mul(V,x,X,a.v,b[i],$$MUL)) + out,t = zeros(svec(N,Any)) .+ out,Any + @inbounds exteraddmulti!(V,out,x,X,derive_mul(V,x,X,a.v,b[i],$$MUL)) + end + else @inbounds outeraddblade!(V,out,x,X,derive_mul(V,x,X,a.v,b[i],$$MUL)) end end - return Chain{V,G+L,t}(out) + return μ ? MultiVector{V}(out) : Chain{V,G+L,t}(out) end end end @generated function contraction(a::Chain{V,G,T},b::Chain{V,L,S}) where {V,G,L,T<:$Field,S<:$Field} - Gndims(V) && (return g_zero(V)) + G+L>ndims(V) && diffvars(V)==0 && (return g_zero(V)) if binomial(ndims(w),G)*binomial(ndims(W),L)<(1<0,mixedmode(W)>0 for i ∈ 1:binomial(ndims(w),G) @inbounds v,iai = :(a[$i]),ia[i] x = CA ? dual(V,iai) : iai for j ∈ 1:binomial(ndims(W),L) X = @inbounds CB ? dual(V,ib[j]) : ib[j] - outeraddblade!_pre(V,out,x,X,derive_pre(V,x,X,v,:(b[$j]),$(QuoteNode(MUL)))) + if μ + exteraddmulti!_pre(V,out,x,X,derive_pre(V,x,X,v,:(b[$j]),$(QuoteNode(MUL)))) + else + outeraddblade!_pre(V,out,x,X,derive_pre(V,x,X,v,:(b[$j]),$(QuoteNode(MUL)))) + end end end - return :(Chain{$V,G+L,promote_type(T,S)}($(Expr(:call,:SVector,out...)))) + return if μ + :(MultiVector{$V}($(Expr(:call,:SVector,out...)))) + else + :(Chain{$V,G+L,promote_type(T,S)}($(Expr(:call,:SVector,out...)))) + end else return quote $(insert_expr((:N,:t,:μ),$(QuoteNode(VEC)))...) ia = indexbasis(ndims(w),G) ib = indexbasis(ndims(W),L) - out = zeros($$VEC(N,G+L,t)) + out = zeros(μ $$VEC(N,t) : $$VEC(N,G+L,t)) CA,CB = mixedmode(w)>0,mixedmode(W)>0 for i ∈ 1:binomial(ndims(w),G) @inbounds v,iai = a[i],ia[i] x = CA ? dual(V,iai) : iai v≠0 && for j ∈ 1:binomial(ndims(W),L) X = @inbounds CB ? dual(V,ib[j]) : ib[j] - if @inbounds outeraddblade!(V,out,x,X,derive_mul(V,x,X,v,b[j],$$MUL))&μ - out,t = zeros(svec(N,G+L,promote_type,Any)) .+ out,Any + if μ + if @inbounds exteraddmulti!(V,out,x,X,derive_mul(V,x,X,v,b[j],$$MUL)) + out,t = zeros(svec(N,promote_type,Any)) .+ out,Any + @inbounds exteraddmulti!(V,out,x,X,derive_mul(V,x,X,v,b[j],$$MUL)) + end + else @inbounds outeraddblade!(V,out,x,X,derive_mul(V,x,X,v,b[j],$$MUL)) end end end - return Chain{V,G+L,t}(out) + return μ ? MultiVector{V}(out) : Chain{V,G+L,t}(out) end end end end