There are significant efficiency gains by implementing Projective Geometric Algebra (PGA) in OpenGL shader hardware as described in Steven De Keninck’s recent paper, “Look, Ma, No Matrices!”
Less significant (because it applies to only a few applications), I’d like to implement a more efficient Julia port of his C++ Projective Geometric Algebra reference implementation by efficiently skipping the calculation of some vector elements that are not needed in some applications. For example, the following code detects the intersection of two line segments, using only the first element of the calculated vector. In this application, calculating the other elements (i.e., 2:8) is a waste.
# mwe.jl - Minimal Working Example
# detect intersecting line segments
# define multivector basis names
# 0 denotes projective dimension (e.g., e0 * e0 = 0)
basis = [ # iField
"1" # 1 scalar
"e0" # 2 grade 1 vectors
"e1" # 3
"e2" # 4
"e01" # 5 grade 2 vectors (bivectors)
"e20" # 6
"e12" # 7
"e012"] # 8 pseudoscalar
# define the basis elements
nField = 2^3 # 3 = 2D + 1 dimensions
e0 = zeros(Float32, nField); e0[2] = 1
e1 = zeros(Float32, nField); e1[3] = 1
e2 = zeros(Float32, nField); e2[4] = 1
e01 = zeros(Float32, nField); e01[5] = 1
e20 = zeros(Float32, nField); e20[6] = 1
e12 = zeros(Float32, nField); e12[7] = 1
e012 = zeros(Float32, nField); e012[8] = 1
# regressive product: vee operator (&, \vee)
function Base.:&(a::Vector{Float32},b::Vector{Float32})::Vector{Float32}
res = similar(a)
res[1]=a[1]*b[8]+a[2]*b[7]+a[3]*b[6]+a[4]*b[5]+a[5]*b[4]+a[6]*b[3]+a[7]*b[2]+a[8]*b[1]
res[2]=a[2]*b[8]+a[8]*b[2]+a[5]*b[6]-a[6]*b[5]
res[3]=a[3]*b[8]+a[8]*b[3]-a[5]*b[7]+a[7]*b[5]
res[4]=a[4]*b[8]+a[8]*b[4]+a[6]*b[7]-a[7]*b[6]
res[5]=a[5]*b[8]+a[8]*b[5]
res[6]=a[6]*b[8]+a[8]*b[6]
res[7]=a[7]*b[8]+a[8]*b[7]
res[8]=a[8]*b[8]
return res
end
function point(x::Number, y::Number)::Vector{Float32}
return x*e20 + y*e01 + e12
end
function tst()
for iTest = 1:2
P0A = point(0,0)
P0B = point(2,0)
P1A = point(1,1)
P1B = iTest == 1 ?
point(0,-1) :
point(-2,-1)
L0 = P0A & P0B # Line 0
A0A = L0 & P1A # oriented area about Line 0
A0B = L0 & P1B # oriented area about Line 0
L1 = P1A & P1B # Line 1
A1A = L1 & P0A # oriented area about Line 1
A1B = L1 & P0B # oriented area about Line 1
println("iTest: $iTest")
println("L0: $L0")
println("A0A: $A0A")
println("A0B: $A0B")
println("L1: $L1")
println("A1A: $A1A")
println("A1B: $A1B")
isIntersection =
(A0A[1] * A0B[1] < 0) &&
(A1A[1] * A1B[1] < 0) ? true : false
println("isIntersection: $isIntersection\n")
end # iTest
end # tst()
julia> include("mwe.jl")
tst (generic function with 1 method)
julia> tst()
iTest: 1
L0: Float32[0.0, 0.0, 0.0, -2.0, 0.0, 0.0, 0.0, 0.0]
A0A: Float32[-2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
A0B: Float32[2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
L1: Float32[0.0, 1.0, -2.0, 1.0, 0.0, 0.0, 0.0, 0.0]
A1A: Float32[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
A1B: Float32[-3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
isIntersection: true
iTest: 2
L0: Float32[0.0, 0.0, 0.0, -2.0, 0.0, 0.0, 0.0, 0.0]
A0A: Float32[-2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
A0B: Float32[2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
L1: Float32[0.0, -1.0, -2.0, 3.0, 0.0, 0.0, 0.0, 0.0]
A1A: Float32[-1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
A1B: Float32[-5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
isIntersection: false
In C or C++, I would reverse the order of vector element calculations and use a switch statement (without break lines) so that a single test would jump to the code to calculate the first n (e.g., n=1 in the line segment intersection code above) elements. However, Julia doesn’t have switch statements and it doesn’t seem that Julia’s macros were intended to “ungenerate” code. Is there an efficient way to skip the calculation of some vector elements?