A fast reference implementation of a pretty projective geometric algebra

General Usage
, , ,
15

Welcome to the club! You might be interested in some of these Julia packages:

I will comment a little more about my own package,

in reference to some of the points you asked about.

Multivector types

Instead of operating directly on Vectors, it’s encouraged to make your own type that’s reflective of the mathematical object you’re representing. My package defines a parametric multivector type (as well as a BasisBlade type):

#                   ┌ metric signature/dimension of algebra
#                   │  ┌ grades represented
#                   │  │ ┌ type of components vector
struct Multivector{Sig,K,S} <: AbstractMultivector{Sig}
	comps::S
end

E.g., Multivector{3,1}(rand(3)) is a 1-vector in ℝ^3, and Multivector{Cl(3,0,1), 0:2:4, MVector{8, Int}} is the type of an even-grade multivector in 3D PGA with Int components stored in a fixed-size MVector. You could even represent a scalar–pseudoscalar combination efficiently as Multivector{8, (0, 8)}([2, 3]) — notice this only has two components (whereas the full multivector has 2^8).

This avoids type piracy, and means all the information you need to know about the dimension, grade, and numeric type of multivector is available at compile-time.

Fast implementations

One possible way of achieving good performance like

your implementation of the PGA 3D geometric product 
function Base.:*(a::Vector{Float32},b::Vector{Float32})
  res = Vector{Float32}(undef,nField)
  res[1]=b[1]*a[1]+b[3]*a[3]+b[4]*a[4]+b[5]*a[5]-b[9]*a[9]-b[10]*a[10]-b[11]*a[11]-b[15]*a[15]
  res[2]=b[2]*a[1]+b[1]*a[2]-b[6]*a[3]-b[7]*a[4]-b[8]*a[5]+b[3]*a[6]+b[4]*a[7]+b[5]*a[8]+b[12]*a[9]+b[13]*a[10]+b[14]*a[11]+b[9]*a[12]+b[10]*a[13]+b[11]*a[14]+b[16]*a[15]-b[15]*a[16]
  res[3]=b[3]*a[1]+b[1]*a[3]-b[9]*a[4]+b[10]*a[5]+b[4]*a[9]-b[5]*a[10]-b[15]*a[11]-b[11]*a[15]
  res[4]=b[4]*a[1]+b[9]*a[3]+b[1]*a[4]-b[11]*a[5]-b[3]*a[9]-b[15]*a[10]+b[5]*a[11]-b[10]*a[15]
  res[5]=b[5]*a[1]-b[10]*a[3]+b[11]*a[4]+b[1]*a[5]-b[15]*a[9]+b[3]*a[10]-b[4]*a[11]-b[9]*a[15]
  res[6]=b[6]*a[1]+b[3]*a[2]-b[2]*a[3]-b[12]*a[4]+b[13]*a[5]+b[1]*a[6]-b[9]*a[7]+b[10]*a[8]+b[7]*a[9]-b[8]*a[10]-b[16]*a[11]-b[4]*a[12]+b[5]*a[13]+b[15]*a[14]-b[14]*a[15]-b[11]*a[16]
  res[7]=b[7]*a[1]+b[4]*a[2]+b[12]*a[3]-b[2]*a[4]-b[14]*a[5]+b[9]*a[6]+b[1]*a[7]-b[11]*a[8]-b[6]*a[9]-b[16]*a[10]+b[8]*a[11]+b[3]*a[12]+b[15]*a[13]-b[5]*a[14]-b[13]*a[15]-b[10]*a[16]
  res[8]=b[8]*a[1]+b[5]*a[2]-b[13]*a[3]+b[14]*a[4]-b[2]*a[5]-b[10]*a[6]+b[11]*a[7]+b[1]*a[8]-b[16]*a[9]+b[6]*a[10]-b[7]*a[11]+b[15]*a[12]-b[3]*a[13]+b[4]*a[14]-b[12]*a[15]-b[9]*a[16]
  res[9]=b[9]*a[1]+b[4]*a[3]-b[3]*a[4]+b[15]*a[5]+b[1]*a[9]+b[11]*a[10]-b[10]*a[11]+b[5]*a[15]
  res[10]=b[10]*a[1]-b[5]*a[3]+b[15]*a[4]+b[3]*a[5]-b[11]*a[9]+b[1]*a[10]+b[9]*a[11]+b[4]*a[15]
  res[11]=b[11]*a[1]+b[15]*a[3]+b[5]*a[4]-b[4]*a[5]+b[10]*a[9]-b[9]*a[10]+b[1]*a[11]+b[3]*a[15]
  res[12]=b[12]*a[1]-b[9]*a[2]+b[7]*a[3]-b[6]*a[4]+b[16]*a[5]-b[4]*a[6]+b[3]*a[7]-b[15]*a[8]-b[2]*a[9]+b[14]*a[10]-b[13]*a[11]+b[1]*a[12]+b[11]*a[13]-b[10]*a[14]+b[8]*a[15]-b[5]*a[16]
  res[13]=b[13]*a[1]-b[10]*a[2]-b[8]*a[3]+b[16]*a[4]+b[6]*a[5]+b[5]*a[6]-b[15]*a[7]-b[3]*a[8]-b[14]*a[9]-b[2]*a[10]+b[12]*a[11]-b[11]*a[12]+b[1]*a[13]+b[9]*a[14]+b[7]*a[15]-b[4]*a[16]
  res[14]=b[14]*a[1]-b[11]*a[2]+b[16]*a[3]+b[8]*a[4]-b[7]*a[5]-b[15]*a[6]-b[5]*a[7]+b[4]*a[8]+b[13]*a[9]-b[12]*a[10]-b[2]*a[11]+b[10]*a[12]-b[9]*a[13]+b[1]*a[14]+b[6]*a[15]-b[3]*a[16]
  res[15]=b[15]*a[1]+b[11]*a[3]+b[10]*a[4]+b[9]*a[5]+b[5]*a[9]+b[4]*a[10]+b[3]*a[11]+b[1]*a[15]
  res[16]=b[16]*a[1]+b[15]*a[2]+b[14]*a[3]+b[13]*a[4]+b[12]*a[5]+b[11]*a[6]+b[10]*a[7]+b[9]*a[8]+b[8]*a[9]+b[7]*a[10]+b[6]*a[11]-b[5]*a[12]-b[4]*a[13]-b[3]*a[14]-b[2]*a[15]+b[1]*a[16]
  return res
end

is to have the compiler compute the resulting expression for you and use that.

I do this (still experimental) my defining a generic inefficient geometric product

function geometric_product(A, B) # PSEUDOCODE
    C = 0
    for a in blades(A), b in blades(B)
        C += geometric_product_blades(a, b)
    end
    C
end

which can be used on multivectors with symbolic components (thanks to Symbolics.jl). Then, the resulting expression can be converted into fully-unrolled code like the function you have above — automatically! (This is achieved with @generated functions.) This efficient code is then reused on similar multivector types. This also means you don’t have to write more code to support other algebras or multivectors of different grade combinations.

Pretty operators

The Julia community definitely loves unicode for purposes like this! We can already define operators like you say. I have

a ∧ b = wedge(a, b)
⋅(a, b) = inner(a, b)
⨼(a, b) = lcontract(a, b)
⨽(a, b) = rcontract(a, b)

but you could also use ∨ for the regressive product and other unicode operators like × for the commutator product, etc.

4 Likes