Rosetta Code

Geometric algebra: Difference between revisions

Page Discussion

Geometric algebra (view source)

Revision as of 11:25, 14 March 2020

16 bytes added ,  4 years ago
Rename Perl 6 -> Raku, alphabetize, minor clean-up
(Rename Perl 6 -> Raku, alphabetize, minor clean-up)
Line 53:
<br><br>
 
=={{header|C++ sharp|C#}}==
{{trans|D}}
<lang cpp>#include <algorithm>
#include <iostream>
#include <random>
#include <vector>
 
double uniform01() {
static std::default_random_engine generator;
static std::uniform_real_distribution<double> distribution(0.0, 1.0);
return distribution(generator);
}
 
int bitCount(int i) {
i -= ((i >> 1) & 0x55555555);
i = (i & 0x33333333) + ((i >> 2) & 0x33333333);
i = (i + (i >> 4)) & 0x0F0F0F0F;
i += (i >> 8);
i += (i >> 16);
return i & 0x0000003F;
}
 
double reorderingSign(int i, int j) {
int k = i >> 1;
int sum = 0;
while (k != 0) {
sum += bitCount(k & j);
k = k >> 1;
}
return ((sum & 1) == 0) ? 1.0 : -1.0;
}
 
struct MyVector {
public:
MyVector(const std::vector<double> &da) : dims(da) {
// empty
}
 
double &operator[](size_t i) {
return dims[i];
}
 
const double &operator[](size_t i) const {
return dims[i];
}
 
MyVector operator+(const MyVector &rhs) const {
std::vector<double> temp(dims);
for (size_t i = 0; i < rhs.dims.size(); ++i) {
temp[i] += rhs[i];
}
return MyVector(temp);
}
 
MyVector operator*(const MyVector &rhs) const {
std::vector<double> temp(dims.size(), 0.0);
for (size_t i = 0; i < dims.size(); i++) {
if (dims[i] != 0.0) {
for (size_t j = 0; j < dims.size(); j++) {
if (rhs[j] != 0.0) {
auto s = reorderingSign(i, j) * dims[i] * rhs[j];
auto k = i ^ j;
temp[k] += s;
}
}
}
}
return MyVector(temp);
}
 
MyVector operator*(double scale) const {
std::vector<double> temp(dims);
std::for_each(temp.begin(), temp.end(), [scale](double a) { return a * scale; });
return MyVector(temp);
}
 
MyVector operator-() const {
return *this * -1.0;
}
 
MyVector dot(const MyVector &rhs) const {
return (*this * rhs + rhs * *this) * 0.5;
}
 
friend std::ostream &operator<<(std::ostream &, const MyVector &);
 
private:
std::vector<double> dims;
};
 
std::ostream &operator<<(std::ostream &os, const MyVector &v) {
auto it = v.dims.cbegin();
auto end = v.dims.cend();
 
os << '[';
if (it != end) {
os << *it;
it = std::next(it);
}
while (it != end) {
os << ", " << *it;
it = std::next(it);
}
return os << ']';
}
 
MyVector e(int n) {
if (n > 4) {
throw new std::runtime_error("n must be less than 5");
}
 
auto result = MyVector(std::vector<double>(32, 0.0));
result[1 << n] = 1.0;
return result;
}
 
MyVector randomVector() {
auto result = MyVector(std::vector<double>(32, 0.0));
for (int i = 0; i < 5; i++) {
result = result + MyVector(std::vector<double>(1, uniform01())) * e(i);
}
return result;
}
 
MyVector randomMultiVector() {
auto result = MyVector(std::vector<double>(32, 0.0));
for (int i = 0; i < 32; i++) {
result[i] = uniform01();
}
return result;
}
 
int main() {
for (int i = 0; i < 5; i++) {
for (int j = 0; j < 5; j++) {
if (i < j) {
if (e(i).dot(e(j))[0] != 0.0) {
std::cout << "Unexpected non-null scalar product.";
return 1;
} else if (i == j) {
if (e(i).dot(e(j))[0] == 0.0) {
std::cout << "Unexpected null scalar product.";
}
}
}
}
}
 
auto a = randomMultiVector();
auto b = randomMultiVector();
auto c = randomMultiVector();
auto x = randomVector();
 
// (ab)c == a(bc)
std::cout << ((a * b) * c) << '\n';
std::cout << (a * (b * c)) << "\n\n";
 
// a(b+c) == ab + ac
std::cout << (a * (b + c)) << '\n';
std::cout << (a * b + a * c) << "\n\n";
 
// (a+b)c == ac + bc
std::cout << ((a + b) * c) << '\n';
std::cout << (a * c + b * c) << "\n\n";
 
// x^2 is real
std::cout << (x * x) << '\n';
 
return 0;
}</lang>
{{out}}
<pre>[-5.63542, -6.59107, -10.2043, -5.21095, 8.68946, 0.579114, -4.67295, -6.72461, 1.55005, -2.63952, -1.83855, 2.4967, -4.3396, -9.9157, -4.6942, -3.23625, -2.3767, -4.55607, -14.3135, -14.2001, 9.84839, 3.69933, -3.38306, -7.60063, -0.236772, 0.988399, -0.549176, 6.61959, 4.69712, -5.34606, -12.2294, -12.6537]
[-5.63542, -6.59107, -10.2043, -5.21095, 8.68946, 0.579114, -4.67295, -6.72461, 1.55005, -2.63952, -1.83855, 2.4967, -4.3396, -9.9157, -4.6942, -3.23625, -2.3767, -4.55607, -14.3135, -14.2001, 9.84839, 3.69933, -3.38306, -7.60063, -0.236772, 0.988399, -0.549176, 6.61959, 4.69712, -5.34606, -12.2294, -12.6537]
 
[-6.27324, -6.7904, 3.10486, -1.14265, 6.38677, 3.57612, -3.90542, -4.17752, -1.36656, -0.0780159, 6.77775, 6.39118, 1.5939, 1.04175, 8.18152, 2.72047, -3.59085, -5.1028, 2.62711, -1.41586, 5.84934, 4.25817, 1.1197, 0.123976, -2.04301, -1.81806, 4.87518, 6.67182, 2.91358, 0.252558, 6.15595, 1.1159]
[-6.27324, -6.7904, 3.10486, -1.14265, 6.38677, 3.57612, -3.90542, -4.17752, -1.36656, -0.0780159, 6.77775, 6.39118, 1.5939, 1.04175, 8.18152, 2.72047, -3.59085, -5.1028, 2.62711, -1.41586, 5.84934, 4.25817, 1.1197, 0.123976, -2.04301, -1.81806, 4.87518, 6.67182, 2.91358, 0.252558, 6.15595, 1.1159]
 
[-7.29133, -8.2555, 0.985588, -1.48171, 2.4995, 4.5152, -1.1938, -2.29702, -2.34025, -2.16526, 10.2208, 7.0629, 0.552639, -0.437582, 7.18962, 2.63274, -3.25348, -4.07006, -0.883786, -3.09677, 1.1018, 2.91198, -0.0095405, 0.0123323, -2.69156, -1.30815, 3.36179, 3.26852, 3.09518, -0.166247, 6.74016, 3.20827]
[-7.29133, -8.2555, 0.985588, -1.48171, 2.4995, 4.5152, -1.1938, -2.29702, -2.34025, -2.16526, 10.2208, 7.0629, 0.552639, -0.437582, 7.18962, 2.63274, -3.25348, -4.07006, -0.883786, -3.09677, 1.1018, 2.91198, -0.0095405, 0.0123323, -2.69156, -1.30815, 3.36179, 3.26852, 3.09518, -0.166247, 6.74016, 3.20827]
 
[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, 0]</pre>
 
=={{header|C#|C sharp}}==
{{trans|D}}
<lang csharp>using System;
Line 429 ⟶ 247:
[-6.85152530876464, -2.6466613868296, -0.814269203555543, -2.63807771948643, 6.89540453623166, 4.70804372948965, 2.78512631505336, 2.03877626337364, -1.77047482836463, 1.52836763170898, 0.71565745873906, 0.886707645932111, -2.25661460785133, -1.08891196061849, 3.44949857952115, 5.8645384965889, -3.09249704978979, -1.41518183096078, 1.87603297737169, 2.45042504250642, 3.66908389117503, 1.85358883025892, 1.23206155683761, -0.216105143607701, -1.88866851071821, -0.131570581491294, 5.7355274883507, 4.22029797577044, -0.212906215521922, -0.323884694190184, 4.41630150883192, 5.94513377054442]
[-6.85152530876464, -2.6466613868296, -0.814269203555542, -2.63807771948643, 6.89540453623166, 4.70804372948965, 2.78512631505335, 2.03877626337364, -1.77047482836463, 1.52836763170898, 0.715657458739061, 0.886707645932111, -2.25661460785133, -1.08891196061849, 3.44949857952115, 5.8645384965889, -3.09249704978979, -1.41518183096078, 1.87603297737169, 2.45042504250642, 3.66908389117503, 1.85358883025892, 1.23206155683761, -0.216105143607701, -1.88866851071821, -0.131570581491294, 5.7355274883507, 4.22029797577044, -0.212906215521922, -0.323884694190183, 4.41630150883192, 5.94513377054442]
 
[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, 0]</pre>
 
=={{header|C++}}==
{{trans|D}}
<lang cpp>#include <algorithm>
#include <iostream>
#include <random>
#include <vector>
 
double uniform01() {
static std::default_random_engine generator;
static std::uniform_real_distribution<double> distribution(0.0, 1.0);
return distribution(generator);
}
 
int bitCount(int i) {
i -= ((i >> 1) & 0x55555555);
i = (i & 0x33333333) + ((i >> 2) & 0x33333333);
i = (i + (i >> 4)) & 0x0F0F0F0F;
i += (i >> 8);
i += (i >> 16);
return i & 0x0000003F;
}
 
double reorderingSign(int i, int j) {
int k = i >> 1;
int sum = 0;
while (k != 0) {
sum += bitCount(k & j);
k = k >> 1;
}
return ((sum & 1) == 0) ? 1.0 : -1.0;
}
 
struct MyVector {
public:
MyVector(const std::vector<double> &da) : dims(da) {
// empty
}
 
double &operator[](size_t i) {
return dims[i];
}
 
const double &operator[](size_t i) const {
return dims[i];
}
 
MyVector operator+(const MyVector &rhs) const {
std::vector<double> temp(dims);
for (size_t i = 0; i < rhs.dims.size(); ++i) {
temp[i] += rhs[i];
}
return MyVector(temp);
}
 
MyVector operator*(const MyVector &rhs) const {
std::vector<double> temp(dims.size(), 0.0);
for (size_t i = 0; i < dims.size(); i++) {
if (dims[i] != 0.0) {
for (size_t j = 0; j < dims.size(); j++) {
if (rhs[j] != 0.0) {
auto s = reorderingSign(i, j) * dims[i] * rhs[j];
auto k = i ^ j;
temp[k] += s;
}
}
}
}
return MyVector(temp);
}
 
MyVector operator*(double scale) const {
std::vector<double> temp(dims);
std::for_each(temp.begin(), temp.end(), [scale](double a) { return a * scale; });
return MyVector(temp);
}
 
MyVector operator-() const {
return *this * -1.0;
}
 
MyVector dot(const MyVector &rhs) const {
return (*this * rhs + rhs * *this) * 0.5;
}
 
friend std::ostream &operator<<(std::ostream &, const MyVector &);
 
private:
std::vector<double> dims;
};
 
std::ostream &operator<<(std::ostream &os, const MyVector &v) {
auto it = v.dims.cbegin();
auto end = v.dims.cend();
 
os << '[';
if (it != end) {
os << *it;
it = std::next(it);
}
while (it != end) {
os << ", " << *it;
it = std::next(it);
}
return os << ']';
}
 
MyVector e(int n) {
if (n > 4) {
throw new std::runtime_error("n must be less than 5");
}
 
auto result = MyVector(std::vector<double>(32, 0.0));
result[1 << n] = 1.0;
return result;
}
 
MyVector randomVector() {
auto result = MyVector(std::vector<double>(32, 0.0));
for (int i = 0; i < 5; i++) {
result = result + MyVector(std::vector<double>(1, uniform01())) * e(i);
}
return result;
}
 
MyVector randomMultiVector() {
auto result = MyVector(std::vector<double>(32, 0.0));
for (int i = 0; i < 32; i++) {
result[i] = uniform01();
}
return result;
}
 
int main() {
for (int i = 0; i < 5; i++) {
for (int j = 0; j < 5; j++) {
if (i < j) {
if (e(i).dot(e(j))[0] != 0.0) {
std::cout << "Unexpected non-null scalar product.";
return 1;
} else if (i == j) {
if (e(i).dot(e(j))[0] == 0.0) {
std::cout << "Unexpected null scalar product.";
}
}
}
}
}
 
auto a = randomMultiVector();
auto b = randomMultiVector();
auto c = randomMultiVector();
auto x = randomVector();
 
// (ab)c == a(bc)
std::cout << ((a * b) * c) << '\n';
std::cout << (a * (b * c)) << "\n\n";
 
// a(b+c) == ab + ac
std::cout << (a * (b + c)) << '\n';
std::cout << (a * b + a * c) << "\n\n";
 
// (a+b)c == ac + bc
std::cout << ((a + b) * c) << '\n';
std::cout << (a * c + b * c) << "\n\n";
 
// x^2 is real
std::cout << (x * x) << '\n';
 
return 0;
}</lang>
{{out}}
<pre>[-5.63542, -6.59107, -10.2043, -5.21095, 8.68946, 0.579114, -4.67295, -6.72461, 1.55005, -2.63952, -1.83855, 2.4967, -4.3396, -9.9157, -4.6942, -3.23625, -2.3767, -4.55607, -14.3135, -14.2001, 9.84839, 3.69933, -3.38306, -7.60063, -0.236772, 0.988399, -0.549176, 6.61959, 4.69712, -5.34606, -12.2294, -12.6537]
[-5.63542, -6.59107, -10.2043, -5.21095, 8.68946, 0.579114, -4.67295, -6.72461, 1.55005, -2.63952, -1.83855, 2.4967, -4.3396, -9.9157, -4.6942, -3.23625, -2.3767, -4.55607, -14.3135, -14.2001, 9.84839, 3.69933, -3.38306, -7.60063, -0.236772, 0.988399, -0.549176, 6.61959, 4.69712, -5.34606, -12.2294, -12.6537]
 
[-6.27324, -6.7904, 3.10486, -1.14265, 6.38677, 3.57612, -3.90542, -4.17752, -1.36656, -0.0780159, 6.77775, 6.39118, 1.5939, 1.04175, 8.18152, 2.72047, -3.59085, -5.1028, 2.62711, -1.41586, 5.84934, 4.25817, 1.1197, 0.123976, -2.04301, -1.81806, 4.87518, 6.67182, 2.91358, 0.252558, 6.15595, 1.1159]
[-6.27324, -6.7904, 3.10486, -1.14265, 6.38677, 3.57612, -3.90542, -4.17752, -1.36656, -0.0780159, 6.77775, 6.39118, 1.5939, 1.04175, 8.18152, 2.72047, -3.59085, -5.1028, 2.62711, -1.41586, 5.84934, 4.25817, 1.1197, 0.123976, -2.04301, -1.81806, 4.87518, 6.67182, 2.91358, 0.252558, 6.15595, 1.1159]
 
[-7.29133, -8.2555, 0.985588, -1.48171, 2.4995, 4.5152, -1.1938, -2.29702, -2.34025, -2.16526, 10.2208, 7.0629, 0.552639, -0.437582, 7.18962, 2.63274, -3.25348, -4.07006, -0.883786, -3.09677, 1.1018, 2.91198, -0.0095405, 0.0123323, -2.69156, -1.30815, 3.36179, 3.26852, 3.09518, -0.166247, 6.74016, 3.20827]
[-7.29133, -8.2555, 0.985588, -1.48171, 2.4995, 4.5152, -1.1938, -2.29702, -2.34025, -2.16526, 10.2208, 7.0629, 0.552639, -0.437582, 7.18962, 2.63274, -3.25348, -4.07006, -0.883786, -3.09677, 1.1018, 2.91198, -0.0095405, 0.0123323, -2.69156, -1.30815, 3.36179, 3.26852, 3.09518, -0.166247, 6.74016, 3.20827]
 
[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, 0]</pre>
Line 1,735:
 
(2.6501903002573224, 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.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)</pre>
 
=={{header|Perl 6}}==
 
Here we write a simplified version of the [https://github.com/grondilu/clifford 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;
}
 
#########################################
## Test code to verify the solution: ##
#########################################
 
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';</lang>
 
=={{header|Phix}}==
Line 2,075 ⟶ 1,970:
 
[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.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, 0.0]</pre>
 
=={{header|Raku}}==
(formerly Perl 6)
 
Here we write a simplified version of the [https://github.com/grondilu/clifford 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;
}
 
#########################################
## Test code to verify the solution: ##
#########################################
 
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';</lang>
 
=={{header|Visual Basic .NET}}==
Thundergnat
10,333

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/262387"