Geometric algebra: Difference between revisions

Content added Content deleted

Inline

Revision as of 22:29, 15 December 2018

Geometric algebra is an other name for Clifford algebras and it's basically an algebra containing a vector space ${\mathcal {V}}$ and obeying the following axioms:

{\begin{array}{c}(ab)c=a(bc)\\a(b+c)=ab+ac\\(a+b)c=ac+bc\\\forall \mathbf {x} \in {\mathcal {V}},\,\mathbf {x} ^{2}\in \mathbb {R} \end{array}}

The product operation in such algebra is called the geometric product. Elements are called multivectors, while multivectors in ${\mathcal {V}}$ are just called vectors.

There are a few simple examples of geometric algebras. A trivial one for instance is simply $\mathbb {R}$ , where ${\mathcal {V}}=\mathbb {R}$ . The complex numbers also form a geometric algebra, where the vector space is the one-dimensional space of all purely imaginary numbers. An other example is the space of quaternions, where the vector space is the three-dimensional space of all linear combinations of $(i,j,k)$ .

The purpose of this task is to implement a geometric algebra with a vector space ${\mathcal {V}}$ of dimension n of at least five, but for extra-credit you can implement a version with n arbitrary large. Using a dimension five is useful as it is the dimension required for the so-called conformal model which will be the subject of a derived task.

To ensure the unicity of the solution (that is, up to some isomorphism), we will also restrict ourselves to the so-called euclidean case, where the square of a non-zero vector is positive:

$\forall \mathbf {x} \in {\mathcal {V}},\,\mathbf {x} \neq 0\implies \mathbf {x} ^{2}>0$ .

You can of course, for extra credit, implement the general case. This would require the definition of a parameter for the signature of the metric.

In order to show that your solution uses a vector space of dimension at least five, you will create a function n -> e(n) such that the vectors e(0), e(1), e(2), e(3), e(4) are linearly independent. To do so you will make them orthonormal with the following scalar product:

$\mathbf {x} \cdot \mathbf {y} =(\mathbf {x} \mathbf {y} +\mathbf {y} \mathbf {x} )/2$

The fact that this so-called inner product defines a scalar product is a consequence of the fourth axiom. To see it one just needs to notice the relation:

$\mathbf {x} \mathbf {y} +\mathbf {y} \mathbf {x} =(\mathbf {x} +\mathbf {y} )^{2}-\mathbf {x} ^{2}-\mathbf {y} ^{2}$

Once you'll have shown that your vector space is at least of dimension five, you will show that the axioms are satisfied. For this purpose you will pick three random multivectors a, b and c, along with a random vector $\mathbf {x}$ .

Producing a random vector is easy. Just use a pseudo-random generation function rand and create a vector:

$\mathrm {randomVector} ()=\sum _{i=0}^{4}\mathrm {rand} ()\mathbf {e} _{i}$

Producing a random multivector is slightly more involved. It is known that when the dimension of ${\mathcal {V}}$ is n, then the dimension of the algebra (seen as a vector space with its natural scalar multiplication) is 2ⁿ. This means that for n=5 there is a basis of 2⁵ = 32 basis multivectors from which any multivector can be written as a linear combination. Create such a basis $m_{0},m_{1},\ldots ,m_{31}$ along with a function producting a random multivector:

$\mathrm {randomMultivector} ()=\sum _{i=0}^{31}\mathrm {rand} ()m_{i}$

To summarize, to solve this task you will:

define the inner product of two vectors : $\mathbf {x} \cdot \mathbf {y} =(\mathbf {xy} +\mathbf {yx} )/2$ .
define the function e
verify the orthonormality $\mathbf {e} _{i}\cdot \mathbf {e} _{j}=\delta _{i,j}$ for i, j in $\{0,1,2,3,4\}$ .
create a function returning a random multivector
create a function returning a random vector
verify the axioms for three rarndom multivectors a, b, c and a random vector x.

Optionally, you will repeat the last step a large number of times, in order to increase confidence in the result.

EchoLisp

We build a CGA based upon a generating quadratic form in R^n. The implementation is general enough, that is ei*ei = +/- 1 , and not restricted to 1. The 5 dimension limit comes from the use of 32 bits numbers to generate all permutations 101... , but this could be improved. The multi-vector multiplication is based on a multiplication table 2^n * 2^n , generated once for all.

<lang scheme> (define e-bits (build-vector 32 (lambda(i) (arithmetic-shift 1 i)))) ;; 1,2,4,.. (define (e-index i) ;; index of ei in native vector (if (zero? i) 0 (arithmetic-shift 1 (1- i))))

(define DIM 0) ;; 2^N (define N 0) (define MultTable null) ;; multiplication table eijk * el.. = exyz.. (define SignTable null) ;; sign of products (define Signature null) ;; input quadratic form

return "eijk"

(define( e-print E sign ) (string-append (cond ((= sign 1) " ") ((= sign -1) "- ") (else "")) (if (zero? E) "1" (for/string ((i N)) #:continue (zero? (bitwise-and E (vector-ref e-bits i))) (string-append "e" (1+ i))))))

returns a string a *e1 + b*e2 + .. z*eijk + ..

(define (multi-print V (x)) (for/string ((i DIM)) (set! x (vector-ref V i)) #:continue (zero? x) (string-append " " (if (> x 0) "+" "") x "*" (e-print i 0))))

generates the multiplication table e_i e__k . * e_j e_l ..==> e_u e_v ... E_I and E_J are sets of indices >=1 , increasing order, represented by a 32 bits number

(define (make-mult-table (verbose #f) (result) (swaps) (ej)) (when verbose (writeln 'N= N 'DIM= DIM 'Q= Signature)) (for* ((E_I (in-range 1 DIM))(E_J (in-range 1 DIM))) (set! result E_I) (set! swaps 0) (for ((j DIM)) ; each bit# in E_J (set! ej (vector-ref e-bits j)) #:continue (zero? (bitwise-and ej E_J))

(for((s (in-range (1- N) j -1))) ;; count swaps (when (!zero? (bitwise-and E_I (vector-ref e-bits s))) (set! swaps (1+ swaps))))

(if (zero? (bitwise-and E_I ej)) ;; e_i * e_j (set! result (bitwise-ior result ej)) (begin ;; else e_i * e_i (set! result (bitwise-xor result ej)) (when (= -1 (vector-ref Signature ej)) (set! swaps (1+ swaps))) ))) ;; j loop

(when verbose (writeln (e-print E_I 0) '* (e-print E_J 0) '= (e-print result (if (even? swaps) 1 -1))))

(matrix-set! MultTable E_I E_J result) (matrix-set! SignTable E_I E_J (if (even? swaps) 1 -1)) ))

multivector operations addition is standard vector addition multiplication a b -> c

(define (multi-mult a b) (define c (make-vector DIM 0)) (for* ((i DIM) (j DIM)) #:continue (zero? (vector-ref a i)) #:continue (zero? (vector-ref b j)) (vector-set! c (array-ref MultTable i j) (+ (* (array-ref SignTable i j) (vector-ref a i) (vector-ref b j)) (vector-ref c (array-ref MultTable i j))))) c)

pretty print a • b or a • b • c

(define ( • a b (c #f)) (multi-print (if c (multi-mult a (multi-mult b c)) (multi-mult a b))))

(Eij i j) -> return multi-vector eiej 0 <= i <= n

(define (Eij i j (coeff 1)) (define Eij (make-vector DIM)) (vector-set! Eij (array-ref MultTable (e-index i) (e-index j)) coeff) Eij)

Reference

https://en.wikipedia.org/wiki/Clifford_algebra#Real_numbers

(make-cga m p [verbose]) => Algebra A(m p) Input a quadratic form Q(x) = x1*x1 + + xm*xm - xm+1*xm+1 - xm+p*xm+p n = m + p = dimension of vector space R^n generates an algebra A(m p) of dimension DIM = 2^n Ex A(n 0) = use R^n dot product as quadratic form : ei*ei = 1 Ex A (0 1) = Complex , e1*e1 = -1 ; A(0 2) => quaternions ei*ei = -1 Implementation limitation n <= 5 multivectors of A(m p) will be mapped on Vectors V of dimension 2^n V[0] is the scalar part of a multivector. Blade of vectors of R^n :V[2^(i-1)] = 1 , 0 elsewhere , i in [1 ..n]

(define (make-cga m p (verbose #f)) (string-delimiter "") (set! N (+ m p)) (set! DIM (expt 2 N)) (set! MultTable (build-array DIM DIM (lambda(i j) (cond ((zero? i) j)((zero? j) i)(else 0))))) (set! SignTable (make-array DIM DIM 1)) (set! Signature (make-vector DIM 1)) ;; Q polynomial (for ((j (in-range m N))) (vector-set! Signature (vector-ref e-bits j) -1))

(make-mult-table verbose) DIM )

</lang>

Output:

we use indices (1 ... n) in conformity with the Wikipedia reference

dimension 2 (1 + e1e2) • (e1 -2e2) = -e1 -3e2

(make-cga 2 0) (define u #(1 0 0 1)) (define v #(0 1 -2 0))

(multi-print u) → +1*1 +1*e1e2 (multi-print v) → +1*e1 -2*e2 (• u v) → -1*e1 -3*e2

task

(make-cga 5 0) (define X #(0 1 -1 0 2 0 0 0 -3 0 0 0 0 0 0 0 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 )) (multi-print X) → +1*e1 -1*e2 +2*e3 -3*e4 -2*e5 (• X X) → +19*1

with another polynomial

(make-cga 0 5) Signature → #( 1 -1 -1 1 -1 1 1 1 -1 1 1 1 1 1 1 1 -1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1) (• X X) → -19*1

(make-cga 4 0) (define i (Eij 1 2)) (define j (Eij 2 3)) (define k (Eij 1 3))

(multi-print i) → +1*e1e2 (multi-print j) → +1*e2e3 (multi-print k)]→ +1*e1e3 (• i i) → -1*1 (• j j) → -1*1 (• k k) → -1*1 (• i j k) → -1*1

(define I (Eij 2 3)) 😖️ error: define : cannot redefine : I (used in Complex) ;; use II instead

(define II (Eij 2 3)) → +1*e2e3 (define J (Eij 3 4)) → +1*e3e4 (define K (Eij 2 4)) → +1*e2e4

(• II II) → -1*1 (• J J) → -1*1 (• K K) → -1*1 (• II J K) → -1*1

</lang> Multiplication table for A(3 0)

N= 3 DIM= 8
e1
e1
e1
e1
e1
e1
e1
e2
e2
e2
e2
e2
e2
e2
e1e2
e1e2
e1e2
e1e2
e1e2
e1e2
e1e2
e3
e3
e3
e3
e3
e3
e3
e1e3
e1e3
e1e3
e1e3
e1e3
e1e3
e1e3
e2e3
e2e3
e2e3
e2e3
e2e3
e2e3
e2e3
e1e2e3
e1e2e3
e1e2e3
e1e2e3
e1e2e3
e1e2e3
e1e2e3 Q= #( 1 1 1 1 1 1 1 1) * e1 = 1 * e2 = e1e2 * e1e2 = e2 * e3 = e1e3 * e1e3 = e3 * e2e3 = e1e2e3 * e1e2e3 = e2e3 * e1 = - e1e2 * e2 = 1 * e1e2 = - e1 * e3 = e2e3 * e1e3 = - e1e2e3 * e2e3 = e3 * e1e2e3 = - e1e3 * e1 = - e2 * e2 = e1 * e1e2 = - 1 * e3 = e1e2e3 * e1e3 = - e2e3 * e2e3 = e1e3 * e1e2e3 = - e3 * e1 = - e1e3 * e2 = - e2e3 * e1e2 = e1e2e3 * e3 = 1 * e1e3 = - e1 * e2e3 = - e2 * e1e2e3 = e1e2 * e1 = - e3 * e2 = - e1e2e3 * e1e2 = e2e3 * e3 = e1 * e1e3 = - 1 * e2e3 = - e1e2 * e1e2e3 = e2 * e1 = e1e2e3 * e2 = - e3 * e1e2 = - e1e3 * e3 = e2 * e1e3 = e1e2 * e2e3 = - 1 * e1e2e3 = - e1 * e1 = e2e3 * e2 = - e1e3 * e1e2 = - e3 * e3 = e1e2 * e1e3 = e2 * e2e3 = - e1 * e1e2e3 = - 1

Go

Translation of: JavaScript

<lang go>package main

import (

   "fmt"
   "math/rand"
   "time"

)

type vector []float64

func e(n uint) vector {

   if n > 4 {
       panic("n must be less than 5")
   }
   result := make(vector, 32)
   result[1<<n] = 1.0
   return result

}

func cdot(a, b vector) vector {

   return mul(vector{0.5}, add(mul(a, b), mul(b, a)))

}

func neg(x vector) vector {

   return mul(vector{-1}, x)

}

func bitCount(i int) int {

   i = i - ((i >> 1) & 0x55555555)
   i = (i & 0x33333333) + ((i >> 2) & 0x33333333)
   i = (i + (i >> 4)) & 0x0F0F0F0F
   i = i + (i >> 8)
   i = i + (i >> 16)
   return i & 0x0000003F

}

func reorderingSign(i, j int) float64 {

   i >>= 1
   sum := 0
   for i != 0 {
       sum += bitCount(i & j)
       i >>= 1
   }
   cond := (sum & 1) == 0
   if cond {
       return 1.0
   }
   return -1.0

}

func add(a, b vector) vector {

   result := make(vector, 32)
   copy(result, a)
   for i, _ := range b {
       result[i] += b[i]
   }
   return result

}

func mul(a, b vector) vector {

   result := make(vector, 32)
   for i, _ := range a {
       if a[i] != 0 {
           for j, _ := range b {
               if b[j] != 0 {
                   s := reorderingSign(i, j) * a[i] * b[j]
                   k := i ^ j
                   result[k] += s
               }
           }
       }
   }
   return result

}

func randomVector() vector {

   result := make(vector, 32)
   for i := uint(0); i < 5; i++ {
       result = add(result, mul(vector{rand.Float64()}, e(i)))
   }
   return result

}

func randomMultiVector() vector {

   result := make(vector, 32)
   for i := 0; i < 32; i++ {
       result[i] = rand.Float64()
   }
   return result

}

func main() {

   rand.Seed(time.Now().UnixNano())
   for i := uint(0); i < 5; i++ {
       for j := uint(0); j < 5; j++ {
           if i < j {
               if cdot(e(i), e(j))[0] != 0 {
                   fmt.Println("Unexpected non-null scalar product.")
                   return
               }
           } else if i == j {
               if cdot(e(i), e(j))[0] == 0 {
                   fmt.Println("Unexpected null scalar product.")
               }
           }
       }
   }

   a := randomMultiVector()
   b := randomMultiVector()
   c := randomMultiVector()
   x := randomVector()

   // (ab)c == a(bc)
   fmt.Println(mul(mul(a, b), c))
   fmt.Println(mul(a, mul(b, c)))

   // a(b + c) == ab + ac
   fmt.Println(mul(a, add(b, c)))
   fmt.Println(add(mul(a, b), mul(a, c)))

   // (a + b)c == ac + bc
   fmt.Println(mul(add(a, b), c))
   fmt.Println(add(mul(a, c), mul(b, c)))

   // x² is real
   fmt.Println(mul(x, x))

}</lang>

Output:

Sample output:

[1.6282881498413662 0.2490818261523896 -0.8936755921478269 -1.555477163901491 6.47589688756284 8.705633181089887 -1.28416798750558 -3.0984267307080446 -14.384954438859133 -13.511137120485879 6.071767421804147 1.099627765550034 -0.40641746354718655 -3.3593459129408076 -2.8089033176352967 -3.003641914720827 4.223517526463662 5.7807271315990185 -4.921271185053852 -5.698203073886508 -0.5449956221104395 3.2199941835007997 -0.4168598688210261 -1.6164380014352773 -13.447900615475964 -11.892642419707807 5.484071302025009 2.781324432176212 5.237445180167182 0.4643791234551212 -7.986755945938485 -2.9272187129576714]
[1.628288149841367 0.24908182615239083 -0.8936755921478254 -1.5554771639014895 6.475896887562836 8.705633181089889 -1.2841679875055805 -3.0984267307080433 -14.384954438859129 -13.51113712048588 6.071767421804145 1.099627765550032 -0.40641746354718755 -3.359345912940806 -2.8089033176352958 -3.003641914720825 4.223517526463663 5.7807271315990185 -4.921271185053855 -5.698203073886506 -0.5449956221104393 3.2199941835008032 -0.4168598688210253 -1.6164380014352775 -13.44790061547596 -11.89264241970781 5.484071302025003 2.781324432176209 5.237445180167183 0.464379123455121 -7.986755945938484 -2.9272187129576706]
[-4.652496095413123 -6.651741805769786 -0.18044192849719706 -1.118756694503706 -1.4545868044605725 0.10199724090664991 0.5018587820915257 2.3004721822960024 -1.813996268087529 0.38357415506855985 7.882236705126414 4.377167004918281 0.338317137066833 1.0631923204534859 5.08861779773926 4.611434580371178 -5.277764644396049 -7.961720991272197 -1.27063408303169 -1.2002120748969933 -0.3251154212726659 2.005000622005424 1.0505371909084391 1.9822320823801767 -2.271503682913346 -0.2403902213900877 6.6269980604812275 4.006018365857085 -0.15074367863748028 -0.48557428903338595 5.291057793190274 2.7751733146879394]
[-4.652496095413124 -6.651741805769786 -0.18044192849719662 -1.1187566945037064 -1.454586804460573 0.10199724090665008 0.5018587820915265 2.3004721822960037 -1.8139962680875295 0.38357415506856163 7.882236705126415 4.377167004918281 0.338317137066834 1.0631923204534859 5.088617797739261 4.611434580371178 -5.277764644396049 -7.9617209912721965 -1.2706340830316898 -1.200212074896994 -0.3251154212726654 2.0050006220054257 1.050537190908438 1.9822320823801762 -2.2715036829133437 -0.24039022139008648 6.626998060481226 4.006018365857084 -0.15074367863747984 -0.4855742890333855 5.291057793190275 2.7751733146879403]
[-4.682894668443622 -7.686899263290272 -0.2500585680601418 -0.8779897639638435 2.579501108403806 2.901924320921563 -0.4542430696469483 1.0374201754917136 -2.588607371991848 -3.229485794697976 7.444924967244714 4.93089687101153 -2.0124785310408284 0.46939350205418884 6.568153651157447 5.710192967946121 -1.41530169954937 -4.536345225213684 0.5240731948596796 1.5123734256564463 1.767786974613905 2.9960861930917018 0.969306657461082 1.036924529835111 -1.456640437477983 -1.3085896429723467 2.9678738261068895 1.3051596528834055 -3.8197616441989037 -1.3030506523824616 4.7818448875243895 4.122121121578954]
[-4.682894668443622 -7.686899263290271 -0.2500585680601421 -0.877989763963844 2.5795011084038055 2.9019243209215633 -0.4542430696469487 1.0374201754917127 -2.5886073719918494 -3.2294857946979736 7.444924967244715 4.93089687101153 -2.0124785310408275 0.46939350205418884 6.568153651157447 5.710192967946121 -1.415301699549369 -4.536345225213685 0.5240731948596802 1.512373425656447 1.7677869746139057 2.9960861930917013 0.9693066574610807 1.036924529835109 -1.456640437477983 -1.3085896429723465 2.9678738261068895 1.305159652883405 -3.819761644198904 -1.3030506523824616 4.781844887524388 4.122121121578957]
[2.1206022437357284 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]

J

Sparse arrays give better performance for this task, than dense arrays, and support relatively arbitrary dimensions, but can be a bit quirky because current implementations of J do not support some features for sparse arrays. We add multivectors x and y using x + y. Also, the first element of one of these multivectors represents the "real valued" or "scalar component" of the multivector.

Implementation:

<lang J>NB. indices are signed machine integers vzero=: 1 $.2^31+IF64*32 odim=. 2^.#vzero

ndx01=:1 :0

 NB. indexed update of numeric rank 1 sparse y
 NB. creating rank 2 sparse result 
 NB. using scalar values from x and scalar inds from m
 NB. where x, m are rank 0 or 1
 NB. (this works around a spurious error in sparse handling)
 n=. #x,.m
 x ((i.n),&.> m)} n#,:y

)

NB. specify that all axes are sparse, for better display clean=: (2;i.@#@$) $. ]

gmul=:4 :0"1

 xj=. ,4$.x
 yj=. ,4$.y
 if. 0= xj *&# yj do. vzero return. end.
 b=. (-##:>./0,xj,yj)&{."1@#:
 xb=. b xj
 yb=. b yj
 rj=. ,#.xb~:"1/yb
 s=. ,_1^ ~:/"1 yb *"1/ 0,.}:"1 ~:/\"1 xb
 vzero (~.rj)}~ rj +//. s*,(xj{x)*/yj{y

)

gdot=: (gmul + gmul~) % 2:

obasis=:1 (2^i.odim)ndx01 vzero e=: {&obasis</lang>

Explanation:

We work with sparse vectors of length 2147483647 on 32 bit machines and of length 9223372036854775807 on 64 bit machines. These are the largest representable vector lengths in current J implementations. J allows negative indices and J uses an index value corresponding to the length of the list to indicate value not found when searching. Thus, current J implementations use signed machine integers for performance and correctness reasons and these are the largest vector lengths we can use in J.

Except, for the purpose of this task, we must pretend that these are "multivectors" instead of vectors - task vectors have a length of log2 the multivector length. So technically speaking, we can only represent 30 element (or less) vectors on 32 bit J and 62 element (or less) vectors on 64 bit J. (Actually, we can represent 31 element vectors on 32 bit J and 63 bit vectors on 64 bit J, but there are hypothetical operations involving the last element which perhaps would be hindered by the fact that the multivector length is not 2147483648 or 9223372036854775808 -- but fortunately, none of this is actually relevant.)

Since these multivectors are silly large, and almost all zeros, we use a sparse representation and only concern ourselves with non-zero values, and their indices.

For addition on multivectors we use J's + (and "multivectors" includes "task vectors").

For "geometric multiplication" (whose relationship to geometry seems obscure, at best) we define something analogous to an inner product. More specifically, we look at the base 2 representation of the argument indices to find the result index and to find whether we negate the product of those two values.

For each pair of values in the left and right argument we use the base 2 representation of the indices to determine what to do with the product of the corresponding two non-zero values:

The result index for that product is the bit-wise exclusive-or of the two indices (the corresponding bit in the result index is 1 when the corresponding bits in the argument indices are different, 0 otherwise).

We negate the product depending on the bits of the left and right argument index. For each 1 bit of the right argument index, we count the number of more significant bits in the left argument index. If this total is odd, we negate the product (but not if it's even). Note that we don't actually have to form a sum to do this - exclusive-or is an adequate replacement for sum if we only need to know whether it's odd or even.

Where we have more than one non-zero value pair contributing to the same product index we sum those values. (Of course, this would apply to value pairs whose result is zero, but since typically there's something approaching 85070591730234615847396907784232501249 of those, we don't actually compute all those zeros...)

Task examples:

<lang J> NB. test arbitrary vector being real (and having the specified result)

  clean gmul~ +/ (e 0 1 2 3 4) gmul 1 _1 2 3 _2 (0 ndx01) vzero

0 │ 19

  NB. required orthogonality
  clean gdot&e&>/~i.4

0 0 0 │ 1 1 1 0 │ 1 2 2 0 │ 1 3 3 0 │ 1

  NB. i j k
  i=: 0 gmul&e 1
  j=: 1 gmul&e 2
  k=: 0 gmul&e 2
   
  i gmul i

0 │ _1

  j gmul j

0 │ _1

  k gmul k

0 │ _1

  i gmul j gmul k

0 │ _1

  NB. I J K
  I=: 1 gmul&e 2
  J=: 2 gmul&e 3
  K=: 1 gmul&e 3
   
  I gmul I

0 │ _1

  J gmul J

0 │ _1

  K gmul K

0 │ _1

  I gmul J gmul K

0 │ _1

K-J

10 │ 1 12 │ _1

  I gmul J+K

10 │ 1 12 │ _1</lang>

Note that sparse arrays display as indices | value for elements which are not the default element (0).

So, for example in the part where we check for orthogonality, we are forming a 5 by 5 by 9223372036854775807 array. The first two dimensions correspond to arguments to e and the final dimension is the multivector dimension. A 0 in the multivector dimension means that that's the "real" or "scalar" part of the multivector. And, since the pair of dimensions for the e arguments whose "dot" product are 1 are identical, we know we have an identity matrix.

JavaScript

<lang javascript>var GA = function () {

   function e(n) {

var result = []; result[1 << n] = 1; return result;

   }
   function cdot(a, b) { return mul([0.5], add(mul(a, b), mul(b, a))) }
   function neg(x) { return mul([-1], x) }
   function bitCount(i) {

// Note that unsigned shifting (>>>) is not required. i = i - ((i >> 1) & 0x55555555); i = (i & 0x33333333) + ((i >> 2) & 0x33333333); i = (i + (i >> 4)) & 0x0F0F0F0F; i = i + (i >> 8); i = i + (i >> 16); return i & 0x0000003F;

   }
   function reorderingSign(a, b) {

a >>= 1; var sum = 0; while (a != 0) { sum += bitCount(a & b); a >>= 1; } return (sum & 1) == 0 ? 1 : -1;

   }
   function add(a, b) {

var result = a.slice(0); for (var i in b) { if (result[i]) { result[i] += b[i]; } else { result[i] = b[i]; } } return result;

   }
   function mul(a, b)
   {

   }
   return {

e : e, cdot : cdot, neg : neg, add : add, mul : mul

};

}();</lang>

And then, from the console:

<lang javascript>var e = GA.e, cdot = GA.cdot;

for (var i = 0; i < 5; i++) {

   for (var j = 0; j < 5; j++) {
       if (i < j) {
           if (cdot(e(i), e(j))[0]) { console.log("unexpected non-nul scalar product"); }
       } else if (i === j) {
           if (!cdot(e(i), e(j))[0]) { console.log("unexpected nul scalar product"); }
       }
   }

}

function randomVector() {

   var result = [];
   for (var i = 0; i < 5; i++) { result = GA.add( result, GA.mul([Math.random()], e(i))); }
   return result;

} function randomMultiVector() {

   var result = [];
   for (var i = 0; i < 32; i++) { result[i] = Math.random(); }
   return result;

}

var a = randomMultiVector(), b = randomMultiVector(), c = randomMultiVector(); var x = randomVector();

// (ab)c == a(bc) console.log(GA.mul(GA.mul(a, b), c)); console.log(GA.mul(a, GA.mul(b, c)));

// a(b + c) == ab + ac console.log(GA.mul(a, GA.add(b, c))); console.log(GA.add(GA.mul(a,b), GA.mul(a, c)));

// (a + b)c == ac + bc console.log(GA.mul(GA.add(a, b), c)); console.log(GA.add(GA.mul(a,c), GA.mul(b, c)));

// x² is real console.log(GA.mul(x, x));</lang>

Output:

[-7.834854130554672, -10.179405417124476, 5.696414143584243, -1.4014556169803851, 12.334288331422336, 11.690738709598888, -0.4279888274147221, 6.226618790084965, -10.904144874917206, -5.46919448234424, -5.647472225071031, -2.9801969751721744, -8.284532508545746, -3.3280413654836494, -2.2182526412098493, 0.4191036292473347, 3.0485450100607103, -0.20619687045226742, 2.1369938048939527, 3.730913391951158, 10.929856967963905, 8.301187183717643, -4.874133827873075, 0.7918650606624789, -8.520661635525103, -7.732342981599732, -6.494750491582618, -2.458749173402162, 3.573788336699224, 2.784339193089742, -1.6479372032388944, -0.35120747879544256]
[-7.83485413055467, -10.179405417124475, 5.696414143584248, -1.4014556169803827, 12.334288331422337, 11.690738709598893, -0.4279888274147213, 6.226618790084964, -10.90414487491721, -5.46919448234424, -5.647472225071032, -2.9801969751721726, -8.284532508545746, -3.3280413654836507, -2.218252641209847, 0.41910362924733874, 3.048545010060707, -0.20619687045226748, 2.136993804893955, 3.7309133919511575, 10.929856967963904, 8.301187183717648, -4.8741338278730755, 0.7918650606624811, -8.520661635525107, -7.732342981599734, -6.494750491582625, -2.45874917340216, 3.5737883366992262, 2.7843391930897443, -1.6479372032388935, -0.351207478795442]
[-4.5157935996060425, -3.9762419076273514, -2.653425845411889, -1.2899302330562412, 6.161562884801266, 3.664812215240675, -0.4471521091019873, 2.39303455739218, -1.6486347268701103, 1.156714478904937, 4.5859158357958965, 6.879356425817299, 1.3341425863947358, 5.641350122882839, 6.378155334673649, 6.466962714879142, -3.645688408496504, -1.9659188980662032, 1.3062519818876646, 1.7973392350972788, 2.4770203476100843, 1.258017836002405, 1.3794942194985413, 3.993871627961031, -3.3620439843097127, -0.4228490927003264, 0.27245046364398495, 3.813642689561589, 2.6785051915908604, 5.409359105713415, 2.9578168177883555, 4.425426168284635]
[-4.515793599606042, -3.976241907627351, -2.653425845411889, -1.2899302330562417, 6.161562884801263, 3.664812215240676, -0.44715210910198766, 2.393034557392179, -1.6486347268701103, 1.156714478904937, 4.585915835795897, 6.8793564258172974, 1.3341425863947352, 5.641350122882839, 6.378155334673649, 6.466962714879143, -3.645688408496502, -1.9659188980662032, 1.3062519818876661, 1.7973392350972783, 2.4770203476100843, 1.258017836002407, 1.379494219498544, 3.99387162796103, -3.3620439843097127, -0.42284909270032545, 0.2724504636439853, 3.8136426895615894, 2.67850519159086, 5.409359105713415, 2.9578168177883555, 4.425426168284636]
[-5.8903316026132755, -6.619647679486295, -1.8140191326116537, -2.519531799741982, 6.604158362571294, 6.352401943423508, 0.9412086471616096, 3.719341486246096, -2.209111542028446, 1.9980997124233557, 5.717878641652222, 7.351597777237362, -2.9037939632499974, 1.497897713658653, 6.811544238648882, 5.861907187665564, -3.2638975880372363, -2.2659714695119115, 1.227221599808634, 0.8343365341022846, 2.72461491531054, 2.728833585944902, 2.226404227376565, 3.888097816250177, 0.35867175462798684, 2.3965356477571302, -1.7151608532791172, 1.403673323043394, -2.1441532262277607, 2.5435142440445646, 2.00110597707534, 1.9825972651495558]
[-5.8903316026132755, -6.6196476794862935, -1.8140191326116533, -2.5195317997419817, 6.604158362571292, 6.352401943423505, 0.9412086471616091, 3.719341486246094, -2.209111542028446, 1.9980997124233555, 5.71787864165222, 7.351597777237364, -2.9037939632499974, 1.4978977136586529, 6.81154423864888, 5.861907187665565, -3.263897588037235, -2.2659714695119124, 1.2272215998086353, 0.8343365341022843, 2.72461491531054, 2.7288335859449018, 2.226404227376565, 3.8880978162501783, 0.3586717546279864, 2.396535647757131, -1.715160853279117, 1.4036733230433938, -2.1441532262277603, 2.543514244044564, 2.0011059770753405, 1.9825972651495565]
[3.193752260485546, 3: 0, 5: 0, 6: 0, 9: 0, 10: 0, 12: 0, 17: 0, 18: 0, 20: 0, 24: 0]

Perl 6

Here we write a simplified version of the Clifford module. It is very general as it is of infinite dimension and also contains an anti-euclidean basis @ē in addition to the euclidean basis @e.

<lang perl6>unit class MultiVector; subset UIntHash of MixHash where .keys.all ~~ UInt; has UIntHash $.blades; method narrow { $!blades.keys.any > 0 ?? self !! ($!blades{0} // 0) }

multi method new(Real $x) returns MultiVector { self.new: (0 => $x).MixHash } multi method new(UIntHash $blades) returns MultiVector { self.new: :$blades }

multi method new(Str $ where /^^e(\d+)$$/) { self.new: (1 +< (2*$0)).MixHash } multi method new(Str $ where /^^ē(\d+)$$/) { self.new: (1 +< (2*$0 + 1)).MixHash }

our @e is export = map { MultiVector.new: "e$_" }, ^Inf; our @ē is export = map { MultiVector.new: "ē$_" }, ^Inf;

my sub order(UInt:D $i is copy, UInt:D $j) {

   (state %){$i}{$j} //= do {

my $n = 0; repeat { $i +>= 1; $n += [+] ($i +& $j).polymod(2 xx *); } until $i == 0; $n +& 1 ?? -1 !! 1;

}

multi infix:<+>(MultiVector $A, MultiVector $B) returns MultiVector is export {

   return MultiVector.new: ($A.blades.pairs, |$B.blades.pairs).MixHash;

} multi infix:<+>(Real $s, MultiVector $B) returns MultiVector is export {

   return MultiVector.new: (0 => $s, |$B.blades.pairs).MixHash;

} multi infix:<+>(MultiVector $A, Real $s) returns MultiVector is export { $s + $A }

multi infix:<*>(MultiVector $, 0) is export { 0 } multi infix:<*>(MultiVector $A, 1) returns MultiVector is export { $A } multi infix:<*>(MultiVector $A, Real $s) returns MultiVector is export {

   MultiVector.new: $A.blades.pairs.map({Pair.new: .key, $s*.value}).MixHash

} multi infix:<*>(MultiVector $A, MultiVector $B) returns MultiVector is export {

   MultiVector.new: do for $A.blades -> $a {

|do for $B.blades -> $b { ($a.key +^ $b.key) => [*] $a.value, $b.value, order($a.key, $b.key), |grep +*, ( |(1, -1) xx * Z* ($a.key +& $b.key).polymod(2 xx *) ) }

   }.MixHash

} multi infix:<**>(MultiVector $ , 0) returns MultiVector is export { MultiVector.new } multi infix:<**>(MultiVector $A, 1) returns MultiVector is export { $A } multi infix:<**>(MultiVector $A, 2) returns MultiVector is export { $A * $A } multi infix:<**>(MultiVector $A, UInt $n where $n %% 2) returns MultiVector is export { ($A ** ($n div 2)) ** 2 } multi infix:<**>(MultiVector $A, UInt $n) returns MultiVector is export { $A * ($A ** ($n div 2)) ** 2 }

multi infix:<*>(Real $s, MultiVector $A) returns MultiVector is export { $A * $s } multi infix:</>(MultiVector $A, Real $s) returns MultiVector is export { $A * (1/$s) } multi prefix:<->(MultiVector $A) returns MultiVector is export { return -1 * $A } multi infix:<->(MultiVector $A, MultiVector $B) returns MultiVector is export { $A + -$B } multi infix:<->(MultiVector $A, Real $s) returns MultiVector is export { $A + -$s } multi infix:<->(Real $s, MultiVector $A) returns MultiVector is export { $s + -$A }

multi infix:<==>(MultiVector $A, MultiVector $B) returns Bool is export { $A - $B == 0 } multi infix:<==>(Real $x, MultiVector $A) returns Bool is export { $A == $x } multi infix:<==>(MultiVector $A, Real $x) returns Bool is export {

   my $narrowed = $A.narrow;
   $narrowed ~~ Real and $narrowed == $x;

}</lang>

And here is the code for verifying the solution:

<lang perl6>use MultiVector; use Test;

plan 29;

sub infix:<cdot>($x, $y) { ($x*$y + $y*$x)/2 }

for ^5 X ^5 -> ($i, $j) {

   my $s = $i == $j ?? 1 !! 0;
   ok @e[$i] cdot @e[$j] == $s, "e$i cdot e$j = $s";

} sub random {

   [+] map {
       MultiVector.new:
       :blades(($_ => rand.round(.01)).MixHash)
   }, ^32;

}

my ($a, $b, $c) = random() xx 3;

ok ($a*$b)*$c == $a*($b*$c), 'associativity'; ok $a*($b + $c) == $a*$b + $a*$c, 'left distributivity'; ok ($a + $b)*$c == $a*$c + $b*$c, 'right distributivity'; my @coeff = (.5 - rand) xx 5; my $v = [+] @coeff Z* @e[^5];

ok ($v**2).narrow ~~ Real, 'contraction';