A fast reference implementation of a pretty projective geometric algebra

General Usage
, , ,
1

As shown at the end of this post, I recently ported pga3d.cpp (the bivector.net C++ reference implementation of projective geometric algebra) to Julia because bivector.net does not currently offer a Julia reference implementation. Before I offer my proposed Julia reference implementation of projective geometric algebra to bivector.net, I am first showing it in this post so my design and implementation can be informally reviewed … because I am relatively new to Julia and I am extremely new to geometric algebra. Also, because some visitors to bivector.net may currently be unfamiliar with Julia, this seems like a good opportunity to make a good first impression of Julia by showcasing its capabilities that seem to be a particularly good fit for implementing projective geometric algebra applications.

DESIGN OVERVIEW AND GOALS
Taking a cue from Elon Musk’s management by rhyming principle “the best part is no part”, the proposed Julia reference implementation defines no structures: the best struct is no struct. Instead, Vector{Float32} is the data type for all the multivectors. Changing the dimension of the space (e.g., 2D, 3D, 4D, 5D) changes the length of the Vector{Float32} and changes some of the overloaded operators but does not change any of the geometric algebra equations.

In his SIBGRAPI 2021 presentation, Steven De Keninck demonstrated the advantage of having geometric algebra equations that stay exactly the same while changing the dimensions of the space. He used ganja.js, his javascript implementation of geometric algebra, to simulate the physics of a 2D wire frame square dangling from a thread, a 3D wire frame cube dangling from a thread, and a 4D wire frame tesseract dangling from a thread … all by changing a single variable defining the number of dimensions. And in his SIGGRAPH 2019 presentation, Steven De Keninck demonstrated how projective geometric algebra simplifies calculating object cross sections. He showed the surprisingly simple projective geometric algebra approach to calculating horizontal slices of a model of a chocolate bunny. (This slicing capability is my motivation for working on this project. More specifically, I want a slicer application implemented in Julia that will generate a 3D printer file that defines the shape of the ratchet mechanism on each of the many display pins in a proposed low cost braille/tactile display.)

A FAST IMPLEMENTATION OF PROJECTIVE GEOMETRIC ALGEBRA
The reason that projective geometric algebra code is surprisingly concise is that the majority of the processing is hidden in overloaded operators. For example, here is the overloaded * operator code that calculates the 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

I attempted a vector multiplication approach to calculating this geometric product but, because it runs 8 times slower, it is currently commented out in the Julia code listed at the end of this post.

According to @time, the proposed Julia implementation at the end of this post runs on my laptop in about 2.2 milliseconds:

julia> @time utest()
point          : 1e032 + 1e123
line           : -1e23
plane          : -3e0 + 2e1 + 1e3
rotor          : 0.7071068 + 0.7071068e12
rotated point  : -0.9999999e013 + 0.9999999e123
rotated line   : 0.9999999e31
rotated plane  : -3e0 + -2e2 + 0.9999999e3
point on plane : 0.2e021 + 1.4e032 + 1e123
point on torus : 0.85e032 + 1e123
toStr() test 1 : -1 + 1e0
toStr() test 2 : 1 + -1e0
  0.002227 seconds (708 allocations: 30.234 KiB)

According to @btime, the proposed Julia implementation at the end of this post runs on my laptop in about 19 microseconds when the output to the screen is disabled:

julia> @btime utest(false)
  19.200 μs (479 allocations: 14.75 KiB)

A PRETTY IMPLEMENTATION OF PROJECTIVE GEOMETRIC ALGEBRA
A standard has emerged on representing geometric algebra operators in programming languages:

Name				      Programming Syntax
geometric product		  a * b
outer product (wedge)	  a ^ b
regressive product (vee)  a & b
inner product (dot)		  a | b
dual				      !a
sandwich product		  a >>> b

Following this standard programming syntax increases the portability of projective geometric algebra equations, even between different programming languages. Unfortunately, this standard programming syntax does not match the standard math syntax. This mismatch between the standard math syntax and the standard programming syntax is an unnecessary mental friction and Julia seems to be in a position to remove it. Perhaps some day Julia will offer a way for programmers to read and write their projective geometric algebra equations directly in math syntax and let Julia convert to and from the programming syntax. For example, perhaps some kind of a block macro could look something like this:

@=pga
   #=
   Inside the pga block macro is a block of projective
   geometry algebra equations using Unicode operators 
   such as \wedge and \vee and \thinspace (for geometric
   product) and \bullet (for dot product) and \asteraccent 
   (for dual) and \tilde (for reverse operator in sandwich).
   =#
end

AN IMPLEMENTATION OF PROJECTIVE GEOMETRIC ALGEBRA

# pga.jl : Projective Geometric Algebra
# This is a port of pga3d.cpp, a C++ reference implementation
# of projective geometric algebra, that is available at
# https://bivector.net/tools.html
#
# The implementation of pga3d.cpp and this port of it, 
# is currently limited to projective geometric algebra in
# 3D space. However, I'm looking forward to extending this 
# implementation to projective geometric algebra in 2D, 4D,
# and possibly 5D spaces. That planned extension is the 
# reason this file is named pga.jl instead of pga3d.jl.
#   
using Printf

# define number of dimensions in space, and
# define number of fields in multivectors
nD = 3
if (nD == 1)
	nField = 0 # TODO: make nonzero after implementation
elseif (nD == 2)
	nField = 0 # TODO: make nonzero after implementation
elseif (nD == 3)
	nField = 16
elseif (nD == 4)
	nField = 0 # TODO: make nonzero after implementation
else
	nField = 0
end

if (nD == 3)
	# define multivector field names
	basis = [		# iField
		"1"		#  1
		"e0"	#  2	basis vectors
		"e1"	#  3
		"e2"	#  4
		"e3"	#  5
		"e01"	#  6	basis bivectors
		"e02"	#  7
		"e03"	#  8
		"e12"	#  9
		"e31"	# 10
		"e23"	# 11
		"e021"	# 12	basis trivectors
		"e013"	# 13
		"e032"	# 14
		"e123"	# 15
		"e0123" # 16	pseudoscalar
	]
	
	# define basis vectors
	e0 = zeros(Float32, nField); e0[2] = 1
	e1 = zeros(Float32, nField); e1[3] = 1
	e2 = zeros(Float32, nField); e2[4] = 1
	e3 = zeros(Float32, nField); e3[5] = 1
	
	# operator overloads dependent upon dimension
	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
	
#	my attempt at speeding up this calculation by vectorizing the products
#	resulted, according to @btime utest(false), a 4.5 times slower unit test
#	function Base.:*(a::Vector{Float32},b::Vector{Float32})
#		res = b .* a[1]
#		
#		tmp = b[[3,1,1,9,10,3,4,5,4,5,15,9,10,11,11,15]] .*
#			  a[[3,2,3,3, 3,2,2,2,3,3, 3,2, 2, 2, 3, 2]]
#		tmp[[5,10,12,13,14]] .*= -1
#		res += tmp
#		
#		tmp = b[[4,6,9,1,11,2,12,13,3,15,5,7,8,16,10,14]] .*
#			  a[[4,3,4,4, 4,3, 3, 3,4, 4,4,3,3, 3, 4, 3]]
#		tmp[[2,3,6,8,9,13]] .*= -1
#		res += tmp
#		
#		tmp = b[[5,7,10,11,1,12,2,14,15,3,4,6,16,8,9,13]] .*
#			  a[[5,4, 5, 5,5, 4,4, 4, 5,5,5,4, 4,4,5, 4]]
#		tmp[[2,4,6,7,11,12]] .*= -1
#		res += tmp
#		
#		tmp = b[[9,8,4,3,15,13,14,2,1,11,10,16,6,7,5,12]] .*
#			  a[[9,5,9,9, 9, 5, 5,5,9, 9, 9, 5,5,5,9, 5]]
#		tmp[[1,2,4,5,7,8,10,14]] .*= -1
#		res += tmp
#		
#		tmp = b[[10,3, 5,15, 3,1,9,10,11, 1, 9,4,5,15, 4,11]] .*
#			  a[[10,6,10,10,10,6,6, 6,10,10,10,6,6, 6,10, 6]]
#		tmp[[1,3,4,8,11,12,14]] .*= -1
#		res += tmp
#		
#		tmp = b[[11,4,15, 5, 4,9,1,11,10, 9, 1,3,15,5, 3,10]] .*
#			  a[[11,7,11,11,11,7,7, 7,11,11,11,7, 7,7,11, 7]]
#		tmp[[1,3,5,6,9,13,14]] .*= -1
#		res += tmp
#		
#		tmp = b[[15,5,11,10, 9,10,11,1, 5, 4, 3,15,3,4, 1,9]] .*
#			  a[[15,8,15,15,15, 8, 8,8,15,15,15, 8,8,8,15,8]]
#		tmp[[1,3,4,5,7,12,13]] .*= -1
#		res += tmp
#		
#		tmp = b[[12,7,6,16,2,14,13,8]] .* a[9]
#		tmp[[3,4,5,6]] .*= -1
#		res[[2,6,7,8,12,13,14,16]] .+= tmp
#		
#		tmp = b[[13,8,16,6,14,2,12,7]] .* a[10]
#		tmp[[2,3,6,7]] .*= -1
#		res[[2,6,7,8,12,13,14,16]] .+= tmp
#		
#		tmp = b[[14,16,8,7,13,12,2,6]] .* a[11]
#		tmp[[2,4,5,7]] .*= -1
#		res[[2,6,7,8,12,13,14,16]] .+= tmp
#		
#		tmp = b[[9,4,3,15,1,11,10,5]] .* a[12]
#		tmp[[2,6,8]] .*= -1
#		res[[2,6,7,8,12,13,14,16]] .+= tmp
#		
#		tmp = b[[10,5,15,3,11,1,9,4]] .* a[13]
#		tmp[[4,7,8]] .*= -1
#		res[[2,6,7,8,12,13,14,16]] .+= tmp
#		
#		tmp = b[[11,15,5,4,10,9,1,3]] .* a[14]
#		tmp[[3,5,8]] .*= -1
#		res[[2,6,7,8,12,13,14,16]] .+= tmp
#		
#		tmp = b[[16,14,13,12,8,7,6,2]] .* a[15]
#		tmp[[2,3,4,8]] .*= -1
#		res[[2,6,7,8,12,13,14,16]] .+= tmp
#		
#		tmp = b[[15,11,10,9,5,4,3,1]] .* a[16]
#		tmp[[1,2,3,4,5,6,7]] .*= -1
#		res[[2,6,7,8,12,13,14,16]] .+= tmp
#		
#		return res
#	end

	function Base.:&(a::Vector{Float32},b::Vector{Float32})
		res = Vector{Float32}(undef,nField)
		res[16]=1*(a[16]*b[16])
		res[15]=-1*(a[15]*-1*b[16]+a[16]*b[15]*-1)
		res[14]=-1*(a[14]*-1*b[16]+a[16]*b[14]*-1)
		res[13]=-1*(a[13]*-1*b[16]+a[16]*b[13]*-1)
		res[12]=-1*(a[12]*-1*b[16]+a[16]*b[12]*-1)
		res[11]=1*(a[11]*b[16]+a[14]*-1*b[15]*-1-a[15]*-1*b[14]*-1+a[16]*b[11])
		res[10]=1*(a[10]*b[16]+a[13]*-1*b[15]*-1-a[15]*-1*b[13]*-1+a[16]*b[10])
		res[9]=1*(a[9]*b[16]+a[12]*-1*b[15]*-1-a[15]*-1*b[12]*-1+a[16]*b[9])
		res[8]=1*(a[8]*b[16]+a[13]*-1*b[14]*-1-a[14]*-1*b[13]*-1+a[16]*b[8])
		res[7]=1*(a[7]*b[16]-a[12]*-1*b[14]*-1+a[14]*-1*b[12]*-1+a[16]*b[7])
		res[6]=1*(a[6]*b[16]+a[12]*-1*b[13]*-1-a[13]*-1*b[12]*-1+a[16]*b[6])
		res[5]=1*(a[5]*b[16]-a[8]*b[15]*-1+a[10]*b[14]*-1-a[11]*b[13]*-1-a[13]*-1*b[11]+a[14]*-1*b[10]-a[15]*-1*b[8]+a[16]*b[5])
		res[4]=1*(a[4]*b[16]-a[7]*b[15]*-1-a[9]*b[14]*-1+a[11]*b[12]*-1+a[12]*-1*b[11]-a[14]*-1*b[9]-a[15]*-1*b[7]+a[16]*b[4])
		res[3]=1*(a[3]*b[16]-a[6]*b[15]*-1+a[9]*b[13]*-1-a[10]*b[12]*-1-a[12]*-1*b[10]+a[13]*-1*b[9]-a[15]*-1*b[6]+a[16]*b[3])
		res[2]=1*(a[2]*b[16]+a[6]*b[14]*-1+a[7]*b[13]*-1+a[8]*b[12]*-1+a[12]*-1*b[8]+a[13]*-1*b[7]+a[14]*-1*b[6]+a[16]*b[2])
		res[1]=1*(a[1]*b[16]+a[2]*b[15]*-1+a[3]*b[14]*-1+a[4]*b[13]*-1+a[5]*b[12]*-1+a[6]*b[11]+a[7]*b[10]+a[8]*b[9]+a[9]*b[8]+a[10]*b[7]+a[11]*b[6]-a[12]*-1*b[5]-a[13]*-1*b[4]-a[14]*-1*b[3]-a[15]*-1*b[2]+a[16]*b[1])
		return res
	end
	
	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[12]*a[4]+b[13]*a[5]+b[1]*a[6]-b[16]*a[11]-b[4]*a[12]+b[5]*a[13]-b[11]*a[16]
		res[7]=b[7]*a[1]+b[12]*a[3]-b[14]*a[5]+b[1]*a[7]-b[16]*a[10]+b[3]*a[12]-b[5]*a[14]-b[10]*a[16]
		res[8]=b[8]*a[1]-b[13]*a[3]+b[14]*a[4]+b[1]*a[8]-b[16]*a[9]-b[3]*a[13]+b[4]*a[14]-b[9]*a[16]
		res[9]=b[9]*a[1]+b[15]*a[5]+b[1]*a[9]+b[5]*a[15]
		res[10]=b[10]*a[1]+b[15]*a[4]+b[1]*a[10]+b[4]*a[15]
		res[11]=b[11]*a[1]+b[15]*a[3]+b[1]*a[11]+b[3]*a[15]
		res[12]=b[12]*a[1]+b[16]*a[5]+b[1]*a[12]-b[5]*a[16]
		res[13]=b[13]*a[1]+b[16]*a[4]+b[1]*a[13]-b[4]*a[16]
		res[14]=b[14]*a[1]+b[16]*a[3]+b[1]*a[14]-b[3]*a[16]
		res[15]=b[15]*a[1]+b[1]*a[15]
		res[16]=b[16]*a[1]+b[1]*a[16]
		return res
	end

	function Base.:^(a::Vector{Float32},b::Vector{Float32})
		res = Vector{Float32}(undef,nField)
		res[1]=b[1]*a[1]
		res[2]=b[2]*a[1]+b[1]*a[2]
		res[3]=b[3]*a[1]+b[1]*a[3]
		res[4]=b[4]*a[1]+b[1]*a[4]
		res[5]=b[5]*a[1]+b[1]*a[5]
		res[6]=b[6]*a[1]+b[3]*a[2]-b[2]*a[3]+b[1]*a[6]
		res[7]=b[7]*a[1]+b[4]*a[2]-b[2]*a[4]+b[1]*a[7]
		res[8]=b[8]*a[1]+b[5]*a[2]-b[2]*a[5]+b[1]*a[8]
		res[9]=b[9]*a[1]+b[4]*a[3]-b[3]*a[4]+b[1]*a[9]
		res[10]=b[10]*a[1]-b[5]*a[3]+b[3]*a[5]+b[1]*a[10]
		res[11]=b[11]*a[1]+b[5]*a[4]-b[4]*a[5]+b[1]*a[11]
		res[12]=b[12]*a[1]-b[9]*a[2]+b[7]*a[3]-b[6]*a[4]-b[4]*a[6]+b[3]*a[7]-b[2]*a[9]+b[1]*a[12]
		res[13]=b[13]*a[1]-b[10]*a[2]-b[8]*a[3]+b[6]*a[5]+b[5]*a[6]-b[3]*a[8]-b[2]*a[10]+b[1]*a[13]
		res[14]=b[14]*a[1]-b[11]*a[2]+b[8]*a[4]-b[7]*a[5]-b[5]*a[7]+b[4]*a[8]-b[2]*a[11]+b[1]*a[14]
		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

	function Base.:~(a::Vector{Float32}) # reverse operator
		res = copy(a)
		res[[6 7 8 9 10 11 12 13 14 15]] .*= -1
		return res
	end

	# basis bivectors (planes)
	function plane(a::Float32, b::Float32, c::Float32, d::Float32) 
		return a*e1 + b*e2 + c*e3 + d*e0
	end
	
	# basis trivectors (points)
	e123 = e1^e2^e3
	e032 = e0^e3^e2
	e013 = e0^e1^e3
	e021 = e0^e2^e1
	function point(x::Float32, y::Float32, z::Float32)
		return e123 + x*e032 + y*e013 + z*e021 
	end
	
	# rotate, translate, conjugate
	function rotor(angle::Float32, line::Vector{Float32})
		return cos(angle/2f0) + sin(angle/2f0)*normalize(line)
	end
	
	function translator(dist::Float32, line::Vector{Float32})
		return 1f0 + dist/2f0 * line
	end
	
	function conjugate(a::Vector{Float32})
		res = copy(a)
		res[[2 3 4 5 6  8 9 10 11]] .*= -1
		return res
	end

	# circle and torus
	function circle(t::Float32, radius::Float32, line::Vector{Float32})
		return rotor(t*2f0*pi,line) * translator(radius,e1*e0)
	end
	
	function torus(s::Float32, t::Float32, 	
		r1::Float32, l1::Vector{Float32},
		r2::Float32, l2::Vector{Float32})
		return circle(s,r2,l2) * circle(t,r1,l1)
	end
end

# overload + operator
function Base.:+(a::Vector{Float32},b::Float32)
	res = copy(a)
	res[1] += b
	return res
end
function Base.:+(a::Float32,b::Vector{Float32})
	res = copy(b)
	res[1] += a
	return res
end

# overload - operator
function Base.:-(a::Vector{Float32},b::Vector{Float32})
	return a - b
end
function Base.:-(a::Vector{Float32},b::Float32)
	res = copy(a)
	res[1] -= b
	return res
end
function Base.:-(a::Float32,b::Vector{Float32})
	res = copy(.-b)
	res[1] += a
	return res
end

# normalization functions
function norm(a::Vector{Float32})
	return sqrt(abs((a * conjugate(a))[1]))
end
function normalize(a::Vector{Float32})
	return a ./ norm(a)
end

# convert multivector fields to string
function toStr(mv::Vector{Float32})
	nField = size(mv,1)
	nNZField = 0
	S = String[]
	for iField = 1:nField
		if mv[iField] != 0
			if nNZField != 0
				push!(S, @sprintf(" + %0.7g%s",
					mv[iField], basis[iField]))
			else
				push!(S, @sprintf("%0.7g%s",
					mv[iField], 
					(iField==1) ? "" : basis[iField]))
			end
			nNZField += 1
		end
	end
	return string(S...)
end

# unit test
# - outputFlag can be set to false to help with benchmarking
# - nD currently does nothing because nD currently is always 3 (like in pga3d.cpp)
function utest(outputFlag::Bool=true, nD::Int64=3)
	axis_z = e1 ^ e2
	origin = axis_z ^ e3
	
	px = point(1f0, 0f0, 0f0)
	if outputFlag
		print("point          : "); println(toStr(px))
	end
	
	line = origin & px
	if outputFlag
		print("line           : "); println(toStr(line))
	end
	
	p = plane(2f0,0f0,1f0,-3f0)
	if outputFlag
		print("plane          : "); println(toStr(p))
	end
	
	rot = rotor(Float32(pi/2), e1*e2)
	if outputFlag
		print("rotor          : "); println(toStr(rot))
	end
	
	rot_point = rot * px * ~rot
	if outputFlag
		print("rotated point  : "); println(toStr(rot_point))
	end
	
	rot_line = rot * line * ~rot
	if outputFlag
		print("rotated line   : "); println(toStr(rot_line))
	end
	
	rot_plane = rot * p * ~rot
	if outputFlag
		print("rotated plane  : "); println(toStr(rot_plane))
	end
	
	point_on_plane = (p | px) * p
	if outputFlag
		print("point on plane : "); println(toStr(normalize(point_on_plane)))
	end
	
	to = torus(0f0, 0f0,
		0.25f0, e1*e2, 
		0.6f0, e1*e3)
	point_on_torus = to * e123 * ~to
	if outputFlag
		print("point on torus : "); println(toStr(point_on_torus))
	end
	
	if outputFlag
		print("toStr() test 1 : "); println(toStr(e0 - 1f0))
		print("toStr() test 2 : "); println(toStr(1f0 - e0))
	end
end
2

Generally, you want to write type-generic code in Jullia. i.e. you should support at least any AbstractVector{<:Real} and probably even AbstractVector{<:Number}. The big advantage of Julia is that you can write type-generic code and it will still be fast, because it is type-specialized at compile-time.

People might want to use other number type double precision, or dual numbers to propagate derivatives.

And people might want to use other vector types. Particular in a small fixed number of dimensions (e.g. 2d, 3d, 4d, or 5d), you almost always want to use StaticArrays.jl for performance reasons.

6 Likes
3

You might have seen some GA in Julia on Github? Perhaps GitHub - chakravala/Grassmann.jl: ⟨Leibniz-Grassmann-Clifford⟩ differential geometric algebra / multivector simplicial complex?

4

I did look at a few Julia geometric algebra packages but my primary interest is in applications and I didn’t see a lot of geometric algebra applications that used a Julia geometric algebra package. In contrast, Steven De Keninck has an impressive collection of geometric algebra applications at https://enkimute.github.io/ganja.js/examples/coffeeshop.html using his ganja.js implementation of geometric algebra. However, perhaps I overlooked something. Are you aware of a Julia geometric algebra package that is being used in a similar collection of freely available geometric algebra applications?

5

I don’t know yet if I’ll encounter an application that requires a lot of multivectors. The projective geometric algebra applications I’ve seen so far are very short.

One could reasonably argue that the projective geometric algebra applications I’ve seen so far are just toy applications meant for demonstrations. However, I think I could productively play with those toy applications for several years. For example, Steven De Keninck demonstrated a one page projective geometric algebra application that calculated the reverse kinematics for a multiple link robot arm to reach some location. I would love to play around with that reverse kinematics application to see if/how projective geometric algebra can handle a severely underconstrained problem (e.g., a robot arm with 20 links).

6

This is type-piracy.
I will let you read the docs on why it is considered poor practice

4 Likes
7

Thanks for the feedback. Is that operator overloading still considered poor practice even when the non-overloaded operation is not allowed? For example, a * b does not compile but a .* b does.

1 Like
8

This case, where you define a method such that it used to be an error and now it does something sensible, is considered a less bad form of piracy; vs change changing the behavior of something that currently works to make it work differently (i.e. monkey patching a special case).

I sometimes call this case, misdemeanor type piracy, and the other felony type piracy.
idk if anyone else uses that terminology.

It is still considered bad, as its kind of overstepping your package’s place.
You are changing someone else’s function on someone elses type.
This can only be done once so if someone else wants to do it, your package and theirs can’t be loaded at the same time.
So as a rule we kinda settle that by saying only the defining package can do this

3 Likes
9

Expecting some sort of warning or error message, I successfully ran the following code in the same Julia REPL session that successfully ran pga.jl, the proposed projective geometric algebra implementation.

# base_conflict.jl : just a test to see how Julia handles
# conflicts in overloading an operator in Base
#   

function Base.:*(a::Vector{Float32},b::Vector{Float32})
	return zeros(Float32,16)
end

# unit test
function utest2(outputFlag::Bool=true, nD::Int64=3)
	res = rand(Float32,16) * rand(Float32,16)
end

Both utest() and utest2() show up in the output of varinfo(). When I repeatedly alternate between running each of them, they both appear to run with their own unique overload of the * operator, as if the operator overload in Base is not a shared resource.

What would I need to do to see the potential penalties for my misdemeanor?

10

Imagine if there were a norm that packages could overload Base.:*(a::Vector{Float32},b::Vector{Float32}) whenever they wanted.

So you had PkgA: *(a::V, b::V) = "a string" and PkgB: *(a::V, b::V) = 1

Then if you did

using PkgA
function foo(a::V, b::V)
    a * b
end
using PkgB
foo()

Using PkgB would totally break PkgA.

So the main problem is that if this were a norm in the Julia community that packages were allowed to do “misdemeanor type piracy”, things would break all the time.

1 Like
11

That’s weird. That is not what i would expect to happen at all.

12

Does that mean I can avoid being a type pirate by not putting pga.jl in a package? Arrr. :grinning:

13

OP did *(a::V, b::V) not import Base; Base.*(a::V, b::V), which is why they are seeing the results they are.

I guess this isn’t type piracy per-se, but it does make code not very readable.

1 Like
14

Better to just make your own types!

Try following this guide and see how far it gets you.

1 Like
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
16

My goal concerning unicode characters in geometric algebra equations is to get the portability of the “standard” programming syntax AND the uncluttered look of the “standard” math syntax (e.g., using spaces instead of asterisks in geometric products). Although I am not yet familiar with VSCode, my hope is that the necessary character substitutions can be done in some sort of VSCode customization.

Thanks. Perhaps one day I will have the math chops to be in such a club, but for now I am just porting code from one language to another in order to get a Julia slicer application for a 3D printer. However, I do have a question for you and/or the club. Are you aware of some capability in the ganja.js environment that has promoted the development of so many projective geometric algebra applications? (I am really impressed by Steven De Keninck’s collection of projective geometric algebra demonstrations at https://enkimute.github.io/ganja.js/examples/coffeeshop.html.)

17

If I were able to edit the title of this post, I would change it to “A fast and pretty reference implementation of projective geometric algebra”, adding the word “reference” to indicate it demonstrates functionality. However, there clearly is an interest in the people who were kind enough to review my design and implementation to add abstraction layers and new type names and packaging. Those additions would probably promote the use of projective geometric algebra in the Julia community, which would be great. So feel free to make those additions but my current focus will remain on the “bare bones” reference implementation, which I think will be enough to implement a Julia slicer application for 3D printing.

Another edit I wish I could make to this post’s title is a rearrangement of the words to “A fast reference implementation of a pretty projective geometric algebra” to more clearly indicate that I’m trying to implement the capability to allow programmers to switch between the portable “standard” programming syntax for geometric algebra equations and the less cluttered (and therefore prettier, in my humble opinion) “standard” math syntax for geometric algebra equations.

Thanks for all the feedback.

18

I updated the title according to your wishes.

1 Like
19

This is an addendum to report and fix a mistake that slowed the execution speed of the unit test by a factor of about 4.2.

Before the fix:

julia> @btime utest(false)
  19.200 μs (479 allocations: 14.75 KiB)

After the fix:

julia> @btime utest(false)
  4.571 μs (77 allocations: 8.70 KiB)

My mistake:

		res = Vector{Float32}(undef,nField)

I used that line of code in several of the overloaded operator implementations to allocate stack memory for a result (“res”) vector. While it was efficient to use undef to avoid unnecessarily zeroing out memory that would soon be set by calculations, it was inefficient to use nField that is defined outside the function. And that slow outside-the-function reference was totally unnecessary given that the length of the result vector is implicitly defined in the function’s vector argument(s).

The fix:
Using length() to calculate the result vector length within the overloaded operator functions would be a reasonable fix. However, I noticed that each of the overloaded operator functions contain at least one component that can be calculated by a vector broadcast operation. Initializing the result vector with the result of that vector broadcast operation not only allocates the result vector memory on the stack, it also slightly reduces the rather lengthy list of multiplies and accumulates in the calculation of the result vector. For example,

#	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 now

	function Base.:*(a::Vector{Float32},b::Vector{Float32})
		res = b .* a[1]
		res[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[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[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[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[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[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[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[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[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[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[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[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[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[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[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[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

Here is the current version of pga.jl including the fix:

# pga.jl : reference implementation of Projective Geometric Algebra
#
# This is a Julia port of pga3d.cpp, bivector.net's C++ 
# reference implementation of projective geometric algebra
# that is available at https://bivector.net/tools.html
#
# The implementation of pga3d.cpp, and this port of it, 
# are currently limited to projective geometric algebra in
# 3D space. However, I'm looking forward to extending this 
# implementation to projective geometric algebra in 2D, 4D,
# and possibly 5D spaces. That planned extension is the 
# reason this file is named pga.jl instead of pga3d.jl.
#   
using Printf

# define number of dimensions in space (e.g., 2D, 3D, 4D)
nD = 3
if (nD == 1)
	nField = 0 # TODO: make nonzero after implementation
elseif (nD == 2)
	nField = 0 # TODO: make nonzero after implementation
	
	# define multivector field names
	basis = [	# iField
		"1"		#  1
		"e0"	#  2	basis vectors
		"e1"	#  3
		"e2"	#  4
		"e01"	#  5	basis bivectors
		"e02"	#  6
		"e12"	#  7
		"e021"	#  8	basis trivectors
		"e0123" #  9	pseudoscalar
	]

elseif (nD == 3)
	nField = 16
	
	# define multivector field names
	basis = [	# iField
		"1"		#  1
		"e0"	#  2	basis vectors
		"e1"	#  3
		"e2"	#  4
		"e3"	#  5
		"e01"	#  6	basis bivectors
		"e02"	#  7
		"e03"	#  8
		"e12"	#  9
		"e31"	# 10
		"e23"	# 11
		"e021"	# 12	basis trivectors
		"e013"	# 13
		"e032"	# 14
		"e123"	# 15
		"e0123" # 16	pseudoscalar
	]
	
	# define basis vectors
	e0 = zeros(Float32, nField); e0[2] = 1
	e1 = zeros(Float32, nField); e1[3] = 1
	e2 = zeros(Float32, nField); e2[4] = 1
	e3 = zeros(Float32, nField); e3[5] = 1
	
	# operator overloads dependent upon dimension
	function Base.:*(a::Vector{Float32},b::Vector{Float32})
		res = b .* a[1]
		res[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[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[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[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[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[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[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[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[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[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[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[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[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[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[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[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
	
#	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
	
#	# NOTE:
#	# According to @btime utest(false, this attempt at speeding up this calculation 
#	# by vectorizing multiplications caused a 4.8 times slower unit test.
#	# Therefore, this vectorized approach is commented out. However, some other
#	# vectorized multiplication may be faster, particularly if implemented in a GPU.
#	function Base.:*(a::Vector{Float32},b::Vector{Float32})
#		res = b .* a[1]
#		
#		tmp = b[[3,1,1,9,10,3,4,5,4,5,15,9,10,11,11,15]] .*
#			  a[[3,2,3,3, 3,2,2,2,3,3, 3,2, 2, 2, 3, 2]]
#		tmp[[5,10,12,13,14]] .*= -1
#		res += tmp
#		
#		tmp = b[[4,6,9,1,11,2,12,13,3,15,5,7,8,16,10,14]] .*
#			  a[[4,3,4,4, 4,3, 3, 3,4, 4,4,3,3, 3, 4, 3]]
#		tmp[[2,3,6,8,9,13]] .*= -1
#		res += tmp
#		
#		tmp = b[[5,7,10,11,1,12,2,14,15,3,4,6,16,8,9,13]] .*
#			  a[[5,4, 5, 5,5, 4,4, 4, 5,5,5,4, 4,4,5, 4]]
#		tmp[[2,4,6,7,11,12]] .*= -1
#		res += tmp
#		
#		tmp = b[[9,8,4,3,15,13,14,2,1,11,10,16,6,7,5,12]] .*
#			  a[[9,5,9,9, 9, 5, 5,5,9, 9, 9, 5,5,5,9, 5]]
#		tmp[[1,2,4,5,7,8,10,14]] .*= -1
#		res += tmp
#		
#		tmp = b[[10,3, 5,15, 3,1,9,10,11, 1, 9,4,5,15, 4,11]] .*
#			  a[[10,6,10,10,10,6,6, 6,10,10,10,6,6, 6,10, 6]]
#		tmp[[1,3,4,8,11,12,14]] .*= -1
#		res += tmp
#		
#		tmp = b[[11,4,15, 5, 4,9,1,11,10, 9, 1,3,15,5, 3,10]] .*
#			  a[[11,7,11,11,11,7,7, 7,11,11,11,7, 7,7,11, 7]]
#		tmp[[1,3,5,6,9,13,14]] .*= -1
#		res += tmp
#		
#		tmp = b[[15,5,11,10, 9,10,11,1, 5, 4, 3,15,3,4, 1,9]] .*
#			  a[[15,8,15,15,15, 8, 8,8,15,15,15, 8,8,8,15,8]]
#		tmp[[1,3,4,5,7,12,13]] .*= -1
#		res += tmp
#		
#		tmp = b[[12,7,6,16,2,14,13,8]] .* a[9]
#		tmp[[3,4,5,6]] .*= -1
#		res[[2,6,7,8,12,13,14,16]] .+= tmp
#		
#		tmp = b[[13,8,16,6,14,2,12,7]] .* a[10]
#		tmp[[2,3,6,7]] .*= -1
#		res[[2,6,7,8,12,13,14,16]] .+= tmp
#		
#		tmp = b[[14,16,8,7,13,12,2,6]] .* a[11]
#		tmp[[2,4,5,7]] .*= -1
#		res[[2,6,7,8,12,13,14,16]] .+= tmp
#		
#		tmp = b[[9,4,3,15,1,11,10,5]] .* a[12]
#		tmp[[2,6,8]] .*= -1
#		res[[2,6,7,8,12,13,14,16]] .+= tmp
#		
#		tmp = b[[10,5,15,3,11,1,9,4]] .* a[13]
#		tmp[[4,7,8]] .*= -1
#		res[[2,6,7,8,12,13,14,16]] .+= tmp
#		
#		tmp = b[[11,15,5,4,10,9,1,3]] .* a[14]
#		tmp[[3,5,8]] .*= -1
#		res[[2,6,7,8,12,13,14,16]] .+= tmp
#		
#		tmp = b[[16,14,13,12,8,7,6,2]] .* a[15]
#		tmp[[2,3,4,8]] .*= -1
#		res[[2,6,7,8,12,13,14,16]] .+= tmp
#		
#		tmp = b[[15,11,10,9,5,4,3,1]] .* a[16]
#		tmp[[1,2,3,4,5,6,7]] .*= -1
#		res[[2,6,7,8,12,13,14,16]] .+= tmp
#		
#		return res
#	end
	
	function Base.:&(a::Vector{Float32},b::Vector{Float32})
		res = a .* b[16]
		res[1] += a[2]*b[15]*-1+a[3]*b[14]*-1+a[4]*b[13]*-1+a[5]*b[12]*-1+a[6]*b[11]+a[7]*b[10]+a[8]*b[9]+a[9]*b[8]+a[10]*b[7]+a[11]*b[6]-a[12]*-1*b[5]-a[13]*-1*b[4]-a[14]*-1*b[3]-a[15]*-1*b[2]+a[16]*b[1]
		res[2] += a[6]*b[14]*-1+a[7]*b[13]*-1+a[8]*b[12]*-1+a[12]*-1*b[8]+a[13]*-1*b[7]+a[14]*-1*b[6]+a[16]*b[2]
		res[3] += -a[6]*b[15]*-1+a[9]*b[13]*-1-a[10]*b[12]*-1-a[12]*-1*b[10]+a[13]*-1*b[9]-a[15]*-1*b[6]+a[16]*b[3]
		res[4] += -a[7]*b[15]*-1-a[9]*b[14]*-1+a[11]*b[12]*-1+a[12]*-1*b[11]-a[14]*-1*b[9]-a[15]*-1*b[7]+a[16]*b[4]
		res[5] += -a[8]*b[15]*-1+a[10]*b[14]*-1-a[11]*b[13]*-1-a[13]*-1*b[11]+a[14]*-1*b[10]-a[15]*-1*b[8]+a[16]*b[5]
		res[6] += a[12]*-1*b[13]*-1-a[13]*-1*b[12]*-1+a[16]*b[6]
		res[7] += -a[12]*-1*b[14]*-1+a[14]*-1*b[12]*-1+a[16]*b[7]
		res[8] += a[13]*-1*b[14]*-1-a[14]*-1*b[13]*-1+a[16]*b[8]
		res[9] += a[12]*-1*b[15]*-1-a[15]*-1*b[12]*-1+a[16]*b[9]
		res[10]+= a[13]*-1*b[15]*-1-a[15]*-1*b[13]*-1+a[16]*b[10]
		res[11]+= a[14]*-1*b[15]*-1-a[15]*-1*b[14]*-1+a[16]*b[11]
		res[12]+= a[16]*b[12]
		res[13]+= a[16]*b[13]
		res[14]+= a[16]*b[14]
		res[15]+= a[16]*b[15]
		return res
	end
	
#	function Base.:&(a::Vector{Float32},b::Vector{Float32})
#		res = Vector{Float32}(undef,nField)
#		res[16]=1*(a[16]*b[16])
#		res[15]=-1*(a[15]*-1*b[16]+a[16]*b[15]*-1)
#		res[14]=-1*(a[14]*-1*b[16]+a[16]*b[14]*-1)
#		res[13]=-1*(a[13]*-1*b[16]+a[16]*b[13]*-1)
#		res[12]=-1*(a[12]*-1*b[16]+a[16]*b[12]*-1)
#		res[11]=1*(a[11]*b[16]+a[14]*-1*b[15]*-1-a[15]*-1*b[14]*-1+a[16]*b[11])
#		res[10]=1*(a[10]*b[16]+a[13]*-1*b[15]*-1-a[15]*-1*b[13]*-1+a[16]*b[10])
#		res[9]=1*(a[9]*b[16]+a[12]*-1*b[15]*-1-a[15]*-1*b[12]*-1+a[16]*b[9])
#		res[8]=1*(a[8]*b[16]+a[13]*-1*b[14]*-1-a[14]*-1*b[13]*-1+a[16]*b[8])
#		res[7]=1*(a[7]*b[16]-a[12]*-1*b[14]*-1+a[14]*-1*b[12]*-1+a[16]*b[7])
#		res[6]=1*(a[6]*b[16]+a[12]*-1*b[13]*-1-a[13]*-1*b[12]*-1+a[16]*b[6])
#		res[5]=1*(a[5]*b[16]-a[8]*b[15]*-1+a[10]*b[14]*-1-a[11]*b[13]*-1-a[13]*-1*b[11]+a[14]*-1*b[10]-a[15]*-1*b[8]+a[16]*b[5])
#		res[4]=1*(a[4]*b[16]-a[7]*b[15]*-1-a[9]*b[14]*-1+a[11]*b[12]*-1+a[12]*-1*b[11]-a[14]*-1*b[9]-a[15]*-1*b[7]+a[16]*b[4])
#		res[3]=1*(a[3]*b[16]-a[6]*b[15]*-1+a[9]*b[13]*-1-a[10]*b[12]*-1-a[12]*-1*b[10]+a[13]*-1*b[9]-a[15]*-1*b[6]+a[16]*b[3])
#		res[2]=1*(a[2]*b[16]+a[6]*b[14]*-1+a[7]*b[13]*-1+a[8]*b[12]*-1+a[12]*-1*b[8]+a[13]*-1*b[7]+a[14]*-1*b[6]+a[16]*b[2])
#		res[1]=1*(a[1]*b[16]+a[2]*b[15]*-1+a[3]*b[14]*-1+a[4]*b[13]*-1+a[5]*b[12]*-1+a[6]*b[11]+a[7]*b[10]+a[8]*b[9]+a[9]*b[8]+a[10]*b[7]+a[11]*b[6]-a[12]*-1*b[5]-a[13]*-1*b[4]-a[14]*-1*b[3]-a[15]*-1*b[2]+a[16]*b[1])
#		return res
#	end
	
	function Base.:|(a::Vector{Float32},b::Vector{Float32})
		res = b .* a[1]
		res[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[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[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[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[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[12]*a[4]+b[13]*a[5]+b[1]*a[6]-b[16]*a[11]-b[4]*a[12]+b[5]*a[13]-b[11]*a[16]
		res[7] += b[12]*a[3]-b[14]*a[5]+b[1]*a[7]-b[16]*a[10]+b[3]*a[12]-b[5]*a[14]-b[10]*a[16]
		res[8] += -b[13]*a[3]+b[14]*a[4]+b[1]*a[8]-b[16]*a[9]-b[3]*a[13]+b[4]*a[14]-b[9]*a[16]
		res[9] += b[15]*a[5]+b[1]*a[9]+b[5]*a[15]
		res[10]+= b[15]*a[4]+b[1]*a[10]+b[4]*a[15]
		res[11]+= b[15]*a[3]+b[1]*a[11]+b[3]*a[15]
		res[12]+= b[16]*a[5]+b[1]*a[12]-b[5]*a[16]
		res[13]+= b[16]*a[4]+b[1]*a[13]-b[4]*a[16]
		res[14]+= b[16]*a[3]+b[1]*a[14]-b[3]*a[16]
		res[15]+= b[1]*a[15]
		res[16]+= b[1]*a[16]
		return res
	end
	
#	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[12]*a[4]+b[13]*a[5]+b[1]*a[6]-b[16]*a[11]-b[4]*a[12]+b[5]*a[13]-b[11]*a[16]
#		res[7]=b[7]*a[1]+b[12]*a[3]-b[14]*a[5]+b[1]*a[7]-b[16]*a[10]+b[3]*a[12]-b[5]*a[14]-b[10]*a[16]
#		res[8]=b[8]*a[1]-b[13]*a[3]+b[14]*a[4]+b[1]*a[8]-b[16]*a[9]-b[3]*a[13]+b[4]*a[14]-b[9]*a[16]
#		res[9]=b[9]*a[1]+b[15]*a[5]+b[1]*a[9]+b[5]*a[15]
#		res[10]=b[10]*a[1]+b[15]*a[4]+b[1]*a[10]+b[4]*a[15]
#		res[11]=b[11]*a[1]+b[15]*a[3]+b[1]*a[11]+b[3]*a[15]
#		res[12]=b[12]*a[1]+b[16]*a[5]+b[1]*a[12]-b[5]*a[16]
#		res[13]=b[13]*a[1]+b[16]*a[4]+b[1]*a[13]-b[4]*a[16]
#		res[14]=b[14]*a[1]+b[16]*a[3]+b[1]*a[14]-b[3]*a[16]
#		res[15]=b[15]*a[1]+b[1]*a[15]
#		res[16]=b[16]*a[1]+b[1]*a[16]
#		return res
#	end
	
	function Base.:^(a::Vector{Float32},b::Vector{Float32})
		res = b .* a[1]
		res[2] += b[1]*a[2]
		res[3] += b[1]*a[3]
		res[4] += b[1]*a[4]
		res[5] += b[1]*a[5]
		res[6] += b[3]*a[2]-b[2]*a[3]+b[1]*a[6]
		res[7] += b[4]*a[2]-b[2]*a[4]+b[1]*a[7]
		res[8] += b[5]*a[2]-b[2]*a[5]+b[1]*a[8]
		res[9] += b[4]*a[3]-b[3]*a[4]+b[1]*a[9]
		res[10]+= -b[5]*a[3]+b[3]*a[5]+b[1]*a[10]
		res[11]+= b[5]*a[4]-b[4]*a[5]+b[1]*a[11]
		res[12]+= -b[9]*a[2]+b[7]*a[3]-b[6]*a[4]-b[4]*a[6]+b[3]*a[7]-b[2]*a[9]+b[1]*a[12]
		res[13]+= -b[10]*a[2]-b[8]*a[3]+b[6]*a[5]+b[5]*a[6]-b[3]*a[8]-b[2]*a[10]+b[1]*a[13]
		res[14]+= -b[11]*a[2]+b[8]*a[4]-b[7]*a[5]-b[5]*a[7]+b[4]*a[8]-b[2]*a[11]+b[1]*a[14]
		res[15]+= 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[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
	
#	function Base.:^(a::Vector{Float32},b::Vector{Float32})
#		res = Vector{Float32}(undef,nField)
#		res[1]=b[1]*a[1]
#		res[2]=b[2]*a[1]+b[1]*a[2]
#		res[3]=b[3]*a[1]+b[1]*a[3]
#		res[4]=b[4]*a[1]+b[1]*a[4]
#		res[5]=b[5]*a[1]+b[1]*a[5]
#		res[6]=b[6]*a[1]+b[3]*a[2]-b[2]*a[3]+b[1]*a[6]
#		res[7]=b[7]*a[1]+b[4]*a[2]-b[2]*a[4]+b[1]*a[7]
#		res[8]=b[8]*a[1]+b[5]*a[2]-b[2]*a[5]+b[1]*a[8]
#		res[9]=b[9]*a[1]+b[4]*a[3]-b[3]*a[4]+b[1]*a[9]
#		res[10]=b[10]*a[1]-b[5]*a[3]+b[3]*a[5]+b[1]*a[10]
#		res[11]=b[11]*a[1]+b[5]*a[4]-b[4]*a[5]+b[1]*a[11]
#		res[12]=b[12]*a[1]-b[9]*a[2]+b[7]*a[3]-b[6]*a[4]-b[4]*a[6]+b[3]*a[7]-b[2]*a[9]+b[1]*a[12]
#		res[13]=b[13]*a[1]-b[10]*a[2]-b[8]*a[3]+b[6]*a[5]+b[5]*a[6]-b[3]*a[8]-b[2]*a[10]+b[1]*a[13]
#		res[14]=b[14]*a[1]-b[11]*a[2]+b[8]*a[4]-b[7]*a[5]-b[5]*a[7]+b[4]*a[8]-b[2]*a[11]+b[1]*a[14]
#		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
	
	function Base.:~(a::Vector{Float32}) # reverse operator
		res = copy(a)
		res[[6 7 8 9 10 11 12 13 14 15]] .*= -1
		return res
	end

	# basis bivectors (planes)
	function plane(a::Float32, b::Float32, c::Float32, d::Float32) 
		return a*e1 + b*e2 + c*e3 + d*e0
	end
	
	# basis trivectors (points)
	e123 = e1^e2^e3
	e032 = e0^e3^e2
	e013 = e0^e1^e3
	e021 = e0^e2^e1
	function point(x::Float32, y::Float32, z::Float32)
		return e123 + x*e032 + y*e013 + z*e021 
	end
	
	# rotate, translate, conjugate
	function rotor(angle::Float32, line::Vector{Float32})
		return cos(angle/2f0) + sin(angle/2f0)*normalize(line)
	end
	
	function translator(dist::Float32, line::Vector{Float32})
		return 1f0 + dist/2f0 * line
	end
	
	function conjugate(a::Vector{Float32})
		res = copy(a)
		res[[2 3 4 5 6  8 9 10 11]] .*= -1
		return res
	end

	# circle and torus
	function circle(t::Float32, radius::Float32, line::Vector{Float32})
		return rotor(t*2f0*pi,line) * translator(radius,e1*e0)
	end
	
	function torus(s::Float32, t::Float32, 	
		r1::Float32, l1::Vector{Float32},
		r2::Float32, l2::Vector{Float32})
		return circle(s,r2,l2) * circle(t,r1,l1)
	end
	
elseif (nD == 4)
	nField = 0 # TODO: make nonzero after implementation
elseif (nD == 5)
	nField = 0 # TODO: make nonzero after implementation
else
	nField = 0
end

if (nD == 3)
end

# overload + operator
function Base.:+(a::Vector{Float32},b::Float32)
	res = copy(a)
	res[1] += b
	return res
end
function Base.:+(a::Float32,b::Vector{Float32})
	res = copy(b)
	res[1] += a
	return res
end

# overload - operator
function Base.:-(a::Vector{Float32},b::Vector{Float32})
	return a - b
end
function Base.:-(a::Vector{Float32},b::Float32)
	res = copy(a)
	res[1] -= b
	return res
end
function Base.:-(a::Float32,b::Vector{Float32})
	res = copy(.-b)
	res[1] += a
	return res
end

# normalization functions
function norm(a::Vector{Float32})
	return sqrt(abs((a * conjugate(a))[1]))
end
function normalize(a::Vector{Float32})
	return a ./ norm(a)
end

# convert multivector fields to string
function toStr(mv::Vector{Float32})
	nField = size(mv,1)
	nNZField = 0
	S = String[]
	for iField = 1:nField
		if mv[iField] != 0
			if nNZField != 0
				push!(S, @sprintf(" + %0.7g%s",
					mv[iField], basis[iField]))
			else
				push!(S, @sprintf("%0.7g%s",
					mv[iField], 
					(iField==1) ? "" : basis[iField]))
			end
			nNZField += 1
		end
	end
	return string(S...)
end

## NOTE:
## According to @btime utest(false), this evaluated version
## of utest() is 98 times slower than the non-evaluated version
## (1.846 ms versus 18.9 us). Therefore, this evaluated version
## of utest() is commented out but not removed ... because it is 
## a reminder to avoid calling eval() when execution speed is an
## issue.
##
## unit test
## - flgOutput should be set to false to help with benchmarking
## - nD currently does nothing because nD currently is always 3 (like in pga3d.cpp)
#function utest(flgOutput::Bool=true, nD::Int64=3)
#	nError = 0
#	axis_z = e1 ^ e2
#	origin = axis_z ^ e3
#	long_exp = quote
#		to = torus(0f0, 0f0, 0.25f0, e1*e2, 0.6f0, e1*e3)
#		point_on_torus = to * e123 * ~to
#	end
#
#	E = [ # Expression array; col 1 is GA equation, col 2 is GA to string conversion
#		:(px = point(1f0,0f0,0f0))      	:(toStr(px));
#		:(line = origin & px)				:(toStr(line));
#		:(p = plane(2f0,0f0,1f0,-3f0))		:(toStr(p));
#		:(rot = rotor(Float32(pi/2),e1*e2))	:(toStr(rot));
#		:(rot_point = rot * px * ~rot)		:(toStr(rot_point));
#		:(rot_line = rot * line * ~rot)		:(toStr(rot_line));
#		:(rot_plane = rot * p * ~rot)		:(toStr(rot_plane));
#		:(point_on_plane = (p | px)*p)		:(toStr(normalize(point_on_plane)));
#		long_exp							:(toStr(point_on_torus));
#		:(tst = e0 - 1f0)					:(toStr(tst));
#		:(tst = -1f0 + e0)					:(toStr(tst));
#	]	
#	S = [ # String array; col 1 is test label, col 2 is expected result
#		" point          : "	"1e032 + 1e123";
#		" line           : "	"-1e23";
#		" plane          : "	"-3e0 + 2e1 + 1e3";
#		" rotor          : "	"0.7071068 + 0.7071068e12";
#		" rotated point  : "	"-0.9999999e013 + 0.9999999e123";
#		" rotated line   : "	"0.9999999e31";
#		" rotated plane  : "	"-3e0 + -2e2 + 0.9999999e3";
#		" point on plane : "	"0.2e021 + 1.4e032 + 1e123";
#		" point on torus : "	"0.85e032 + 1e123";
#		" toStr() test 1 : "	"-1 + 1e0";
#		" toStr() test 2 : "	"1 + -1e0";
#	]
#	
#	nTest = size(S,1)
#	for iTest = 1:nTest
#		eval(E[iTest,1])
#		str = eval(E[iTest,2])
#		if flgOutput
#			isError = str != S[iTest,2]
#			xChar = isError ? 'x' : ' '
#			println(xChar * S[iTest,1] * str)
#			if (isError)
#				println(' ' * S[iTest,1] * str)
#				nError += 1
#			end
#		end
#	end
#	return nError # 'x' is first character in lines with errors
#end

# unit test
# - flgOutput should be set to false to help with benchmarking
# - nD currently does nothing because nD currently is always 3 (like in pga3d.cpp)
function utest(flgOutput::Bool=true, nD::Int64=3)
	nError = 0
	
	# geometric algebra equations in programming syntax
	axis_z = e1 ^ e2
	origin = axis_z ^ e3
	
	px = point(1f0, 0f0, 0f0)
	line = origin & px
	p = plane(2f0,0f0,1f0,-3f0)
	rot = rotor(Float32(pi/2), e1*e2)
	rot_point = rot * px * ~rot
	rot_line = rot * line * ~rot
	rot_plane = rot * p * ~rot
	point_on_plane = (p | px) * p
	to = torus(0f0, 0f0,
		0.25f0, e1*e2, 
		0.6f0, e1*e3)
	point_on_torus = to * e123 * ~to
	
	tst1 = e0 - 1f0
	tst2 = 1f0 - e0
	
	# geometric algebra equations in math syntax
	# axis_z = e1 ∧ e2
	
	# if (slow) output of unit test results wanted
	if flgOutput
		S = Matrix{String}(undef,11,3) # 3 columns:
		S[1,1] = " point          : "  # 1) label
		S[1,2] = toStr(px)             # 2) toStr()
		S[1,3] = "1e032 + 1e123"       # 3) expected string
		
		S[2,1] = " line           : "
		S[2,2] = toStr(line)
		S[2,3] = "-1e23"
		
		S[3,1] = " plane          : "
		S[3,2] = toStr(p)
		S[3,3] = "-3e0 + 2e1 + 1e3"
		
		S[4,1] = " rotor          : "
		S[4,2] = toStr(rot)
		S[4,3] = "0.7071068 + 0.7071068e12"
		
		S[5,1] = " rotated point  : "
		S[5,2] = toStr(rot_point)
		S[5,3] = "-0.9999999e013 + 0.9999999e123"
		
		S[6,1] = " rotated line   : "
		S[6,2] = toStr(rot_line)
		S[6,3] = "0.9999999e31"
		
		S[7,1] = " rotated plane  : "
		S[7,2] = toStr(rot_plane)
		S[7,3] = "-3e0 + -2e2 + 0.9999999e3"
		
		S[8,1] = " point on plane : "
		S[8,2] = toStr(normalize(point_on_plane))
		S[8,3] = "0.2e021 + 1.4e032 + 1e123"
		
		S[9,1] = " point on torus : "
		S[9,2] = toStr(point_on_torus)
		S[9,3] = "0.85e032 + 1e123"
		
		S[10,1] = " toStr() test 1 : "
		S[10,2] = toStr(tst1)
		S[10,3] = "-1 + 1e0"
		
		S[11,1] = " toStr() test 2 : "
		S[11,2] = toStr(tst2)
		S[11,3] = "1 + -1e0"
		
		# print unit test results;
		# 'x' in first column in tests with errors
		nTest = size(S,1)
		for iTest = 1:nTest
			isError = S[iTest,2] != S[iTest,3]
			xChar = isError ? 'x' : ' '
			println(xChar * S[iTest,1] * S[iTest,2])
			if (isError)
				println(' ' * S[iTest,1] * S[iTest,3])
				nError += 1
			end
		end
	end
	return nError # return unit test results
end