Geometric algebra

Geometric algebra is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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 this algebra is called geometric product. The elements are called multivectors. Multivectors in ${\mathcal {V}}$ are just called vectors.

The purpose of this task is to implement such an algebra with vectors of arbitrary size, or up to 32 dimensions if that's easier to implement in your language.

To demonstrate your solution, you will use it to implement quaternions:

• define a scalar product of two vectors as the symmetric part of the geometric product $\mathbf {a} \cdot \mathbf {b} =(\mathbf {ab} +\mathbf {ba} )/2$
• define a function $n\mapsto e(n)$ which allows for creation of an orthonormal basis of ${\mathcal {V}}$

$(\mathbf {e} _{0},\mathbf {e} _{1},\mathbf {e} _{2},\ldots )$,

• verify that the square of a random vector is real. Use this example: $(\mathbf {e} _{0}-\mathbf {e} _{1}+2\mathbf {e} _{2}+3\mathbf {e} _{3}-2\mathbf {e} _{4})^{2}=19$
• verify the orthonormality $\mathbf {e} _{i}\cdot \mathbf {e} _{j}=\delta _{i,j}$ for i, j in $\{0,1,2,3\}$.
• define the following three constants:

${\begin{array}{c}i=\mathbf {e} _{0}\mathbf {e} _{1}\\j=\mathbf {e} _{1}\mathbf {e} _{2}\\k=\mathbf {e} _{0}\mathbf {e} _{2}\end{array}}$

• verify that $i^{2}=j^{2}=k^{2}=ijk=-1$
• to verify that your implementation does not just reuse an ad-hoc implementation of quaternions, also define the following three constants:

${\begin{array}{c}I=\mathbf {e} _{1}\mathbf {e} _{2}\\J=\mathbf {e} _{2}\mathbf {e} _{3}\\K=\mathbf {e} _{1}\mathbf {e} _{3}\end{array}}$

• verify that $I^{2}=J^{2}=K^{2}=IJK=-1$

J

Based on the quaternion implementation at quaternion type, but implementing a clifford algebra which specifically corresponds to the task requirements.

Implementation:

<lang J>dim=.32

T=.3 :0 dim

 s=: _1^~:/@,@(* 1 |:@}._1&(|.!.0)^:a:)/ ::0:@#:@,"0/~x: i. y
 s (<@([ , #.@:(~:/)@#:@, , ])"0/~i.y) } (3#y)$0

)

domain=: 32&{. :.(#~ +./\.@:(0&~:))

ip=: +/ .* gmul=: (ip T ip ])&.domain"1 gadd=: +&.domain"1 gdot=: (gmul + gmul~)&.domain % 2:

e=: {&((2^i.2^.dim){=i. dim)</lang>

Use:

<lang J> NB. test arbitrary vector being real

  gmul~ gadd/ (,.1 _1 2 3 _2) gmul e 0 1 2 3 4

19

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

┌─┬─┬─┬─┐ │1│ │ │ │ ├─┼─┼─┼─┤ │ │1│ │ │ ├─┼─┼─┼─┤ │ │ │1│ │ ├─┼─┼─┼─┤ │ │ │ │1│ └─┴─┴─┴─┘

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

_1

  j gmul j

_1

  k gmul k

_1

  i gmul j gmul k

_1

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

_1

  J gmul J

_1

  K gmul K

_1

  I gmul J gmul K

_1</lang>

JavaScript

<lang javascript>var CGA = function () {

   function e(n) {

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

   }
   function neg(x) { return multiply([-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 = []; for (var i = 0; i < 32; i++) { if (a[i] && b[i]) { var r = a[i] + b[i]; if (r !== 0) { result[i] = r; } } else if (a[i]) { result[i] = a[i]; } else if (b[i]) { result[i] = b[i]; } } return result;

   }
   function multiply(a, b)
   {

var result = []; for (var i = 0; i < 32; i++) { if (a[i]) { for (var j = 0; j < 32; j++) { if (b[j]) { var s = reorderingSign(i, j) * a[i] * b[j]; // if (i == 1 && j == 1) { s *= -1 } // e0*e0 == -1 var k = i ^ j; if (result[k]) { result[k] += s; } else { result[k] = s; } } } } } return result;

   }
   return {

e  : e, neg : neg, add : add, mul : multiply

   };

}(); </lang>

And then, from the console:

<lang javascript>var e = CGA.e; function cdot(a, b) { return CGA.mul([0.5], CGA.add(CGA.mul(a, b), CGA.mul(b, a))) }

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

   for (var j = 0; j < 4; 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"); }
       }
   }

}

var v = e(0); v = CGA.add(v, CGA.mul([-1], e(1))); v = CGA.add(v, CGA.mul([2], e(2))); v = CGA.add(v, CGA.mul([3], e(3))); v = CGA.add(v, CGA.mul([-2], e(4)));

console.log(CGA.mul(v, v)); // [19, 3: 0, 5: 0, 6: 0, 9: 0, 10: 0, 12: 0, 17: 0, 18: 0, 20: 0, 24: 0]

var i = CGA.mul(e(0), e(1)); var j = CGA.mul(e(1), e(2)); var k = CGA.mul(e(0), e(2));

console.log(CGA.mul(i, i)); // [-1] console.log(CGA.mul(j, j)); // [-1] console.log(CGA.mul(k, k)); // [-1] console.log(CGA.mul(CGA.mul(i, j), k)); // [-1]

var I = CGA.mul(e(1), e(2)); var J = CGA.mul(e(2), e(3)); var K = CGA.mul(e(1), e(3));

console.log(CGA.mul(I, I)); // [-1] console.log(CGA.mul(J, J)); // [-1] console.log(CGA.mul(K, K)); // [-1] console.log(CGA.mul(CGA.mul(I, J), K)); // [-1]</lang>

Perl 6

<lang perl6>unit class MultiVector; has Real %.blades{UInt}; method clean { for %!blades { %!blades{.key} :delete unless .value; } } method narrow {

   for %!blades { return self if .key > 0 && .value !== 0; }
   return %!blades{0} // 0;

}

sub e(UInt $n?) returns MultiVector is export {

   $n.defined ?? MultiVector.new(:blades(my Real %{UInt} = (1 +< $n) => 1)) !! MultiVector.new

}

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

   my $n = 0;
   repeat {

$i +>= 1; $n += [+] ($i +& $j).base(2).comb;

   } until $i == 0;
   return $n +& 1 ?? -1 !! 1;

}

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

   my Real %blades{UInt} = $A.blades.clone;
   for $B.blades {

%blades{.key} += .value; %blades{.key} :delete unless %blades{.key};

   }
   return MultiVector.new: :%blades;

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

   my Real %blades{UInt} = $A.blades.clone;
   %blades{0} += $s;
   %blades{0} :delete unless %blades{0};
   return MultiVector.new: :%blades;

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

   my Real %blades{UInt};
   for $A.blades -> $a {

for $B.blades -> $b { my $c = $a.key +^ $b.key; %blades{$c} += $a.value * $b.value * order($a.key, $b.key); %blades{$c} :delete unless %blades{$c}; }

   }
   return MultiVector.new: :%blades;

} 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:<*>(MultiVector $, 0) returns MultiVector is export { MultiVector.new } multi infix:<*>(MultiVector $A, 1) returns MultiVector is export { $A } multi infix:<*>(MultiVector $A, Real $s) returns MultiVector is export {

   return MultiVector.new: :blades(my Real %{UInt} = map { .key => $s * .value }, $A.blades);

} 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 implementing and verifying quaternions:

<lang perl6>use MultiVector; use Test;

plan 13;

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

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

   next if $i > $j;
   ok e($i) cdot e($j) == +($i == $j);

}

my $v = e(0) - e(1) + 2*e(2) + 3*e(3) - 2*e(4);

ok $v**2 == 19, 'v² = 19';

my constant i = e(0)*e(1); my constant j = e(1)*e(2); my constant k = e(0)*e(2);

ok i**2 == j**2 == k**2 == i*j*k == -1, "i² = j² = k² = ijk = -1";

my constant I = e(1)*e(2); my constant J = e(2)*e(3); my constant K = e(1)*e(3);

ok I**2 == J**2 == K**2 == I*J*K == -1, "I² = J² = K² = IJK = -1";</lang>

Output:

1..13
ok 1 - 
ok 2 - 
ok 3 - 
ok 4 - 
ok 5 - 
ok 6 - 
ok 7 - 
ok 8 - 
ok 9 - 
ok 10 - 
ok 11 - v² = 19
ok 12 - i² = j² = k² = ijk = -1
ok 13 - I² = J² = K² = IJK = -1