Friday 11 March 2011

A Few Hypercomplex Onions

C# Generic Limitations

The C# class Onion (listing below) allows experimentation with any of the Cayley-Dickson algebras mentioned in the previous related article, from the monodimensional reals, right up to the 512-dimensional austons :-) and beyond.

Modern professional libraries for this purpose still tend to be based on technologies like C++ templates, allowing for type variation in the underlying coefficients. And not just between a few integer and floating point types. For example, with careful design you can define your Biquaternion type as Quaternion<Complex>. But while elegant and fun to develop, these techniques don't port well to C#. There are several reasons for this, two of them being the comparitively restrictive natures of C# generics, and of C# operator overloading.

Introducing The Onion

The obvious alternative is to create a single class to represent a Cayley-Dickson number of any dimension. That's what Onion does. While it has a private array of double called Values, containing its real coefficients, it also makes essential use of the composition build model. Complex numbers can be treated as ordered pairs of real numbers; quaternions as ordered pairs of complex; octonions, pairs of quaternions; and so on. To support this model, an Onion has read only properties named A and B returning these two lower dimensional Onions.

These let the class encapsulate the Cayley-Dickson definition of a product. Take any two Onions, p and q, of the same dimensionality:
p = (a, b)
q = (c, d)
The Cayley-Dickson construction defines the product of these as
pq = (ac - d'b, da + bc')
where c' and d' represent the conjugates of c and d. Since all of a, b, c, d have lower dimension than p and q, we can apply this definition recursively until we reach the reals, which we already knew how to multiply together. The implementation of operator * is complicated by the need to handle input parameters p, q of different dimensionalities; in such cases, the "missing" component (b or d) is replaced by zero.

Onion Operators, Properties, Methods

Operators implemented in Onion include the unary +, -, and also ~, which returns the conjugate; the binary ==, !=, +, -, *, /, and also ^, which has an optimized override raising an Onion to a signed integral power. Finally, there is an implicit conversion operator from type double, which allows real values to be intermixed with Onions in arithmetic expressions.

Care and parentheses are recommended with ^. In normal C# usage, ^ represents the logical XOR operation, and so has an inappropriately low precedence for the present context (midway between logical OR and logical AND). Also, the exponentiation operator is frequently expected to be right-associative:
2 ^ 3 ^ 4 = 2 ^ (3 ^ 4) = 2 ^ 81 = 2,417,851,639,229,258,349,412,352
whereas without explicit parentheses, we get this instead:
2 ^ 3 ^ 4 = (2 ^ 3) ^ 4 = 8 ^ 4 = 4,096
Most other properties and methods are self-explanatory. Dim is the dimension of the Onion; Re and Im its real and imaginary parts; Norm its norm or modulus, Norm_2 the square of that; and Reciprocal - well, take a wild guess. There are static Exp and (natural) Log functions, added for the particular investigations I happened to be interested in at the time; and a useful ToString() override using 1 + 2i + 3j + 4ij notation, omitting redundant 0s and 1s. Note the use of ij in place of the conventional quaternion k here, as in my earlier article; if you want to change the basis representation, private method AppendDim is the place to visit.

It's a Reference Type but...

Onion isn't quite as dynamic as it appears. Every algebra that it can represent is closed. For example, starting with a bunch of complex numbers, there's no way that any sequence of additions, multiplications etc. can generate a quaternion. Also, I've made it immutable, because it's not a struct (does not have value copy semantics) and I didn't want multireferenced instances changing in unexpected ways.

There is an indexer, but that just gives read access to the coefficients. Once created, an Onion cannot change its value. Expressions such as z *= 2, rather than operating on z itself, replace it with a new Onion.

Examples

Here's some sample code showing verification of the "odd" results used in the previous article:
var i = new Onion(0, 1);
Onion p = i, z1 = (p ^ 2) + 1; // Caution! Operator ^ has a low precedence!
Console.WriteLine("(Complex numbers) When p = {0}, then p² + 1 = {1}.", p, z1);
// Output: (Complex numbers) When p = i, then p² + 1 = 0.

var j = new Onion(0, 0, 1);
Onion q = j, z2 = (p*q - q*p) ^ 2;
Console.WriteLine("(Quaternions) When p = {0} and q = {1}, then (pq - qp)² = {2}.", p, q, z2);
// Output: (Quaternions) When p = i and q = j, then (pq - qp)² = -4.

var k = new Onion(0, 0, 0, 0, 1);
Onion r = k, z3 = ((p*q)*r - p*(q*r)) ^ 2;
Console.WriteLine("(Octonions) When p = {0}, q = {1} and r = {2}, then ((pq)r - p(qr))² = {3}.", p, q, r, z3);
// Output: (Octonions) When p = i, q = j and r = k, then ((pq)r - p(qr))² = -4.

var l = new Onion(0, 0, 0, 0, 0, 0, 0, 0, 1);
p = i*j + k*l; // p ≠ 0
q = i*k + j*l; // q ≠ 0
var pq = p*q;
Console.WriteLine("(Sedenions) When p = {0} and q = {1}, then pq = {2};", p, q, pq);
// Output: (Sedenions) When p = ij + kl and q = ik + jl, then pq = 0;
Console.WriteLine("(p == 0) = {0}; (q == 0) = {1}; (pq == 0) = {2}.", p == 0, q == 0, pq == 0);
// Output: (p == 0) = False; (q == 0) = False; (pq == 0) = True.
Another example: here is the code fragment I used to generate the sedenion multiplication table (order = 4, dim = 16) in the previous article.
const int order = 4;
const int dim = 1 << order;
var values = new double[dim];
for (var row = 0; row < dim; row++)
{
values[row] = 1;
var p = new Onion(values);
values[row] = 0;
for (var col = 0; col < dim; col++)
{
values[col] = 1;
var q = new Onion(values);
values[col] = 0;
Console.Write(string.Format("{{0,{0}}}", order + 2), p*q);
}
Console.WriteLine();
}
Source Code

Okay, finally here's class Onion:
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Onions
{
public class Onion
{
private double[] Values { get; set; }

public Onion(Onion x) : this(x.Values) { }

public Onion(Onion a, Onion b)
{
var dim = Math.Max(a.Dim, b.Dim);
Dim = dim*2;
a.Values.CopyTo(Values, 0);
b.Values.CopyTo(Values, dim);
}

public Onion(params double[] values)
{
var count = values.Count();
var dim = 1;
while (dim < count) dim <<= 1;
Dim = dim;
values.CopyTo(Values, 0);
}

private Onion(IEnumerable<double> values) : this(values.ToArray()) { }

public double this[int index] { get { return index < Dim ? Values[index] : 0; } }

public int Dim
{
get { return Values == null ? 0 : Values.Length; }
private set { if (Dim != value) Values = new double[value]; }
}

public Onion A { get { return new Onion(Values.Take(Dim/2)); } }

public Onion B { get { var dim = Dim/2; return new Onion(Values.Skip(dim).Take(dim)); } }

public double Re { get { return this[0]; } }

public Onion Im { get { return this - Re; } }

public double Norm { get { return Dim == 1 ? Math.Abs(Re) : Math.Sqrt(Norm_2); } }

private double Norm_2 { get { return Values.Sum(a => a*a); } }

public Onion Reciprocal { get { return Dim == 1 ? 1/this[0] : ~this/Norm_2; } }

public static implicit operator Onion(double value) { return new Onion(value); }

public static Onion operator +(Onion x) { return x; }

public static Onion operator -(Onion x) { return Operate(x, a => -a); }

public static Onion operator ~(Onion x) { return x.Dim == 1 ? x : new Onion(~x.A, -x.B); }

public static bool operator ==(Onion x, Onion y)
{
bool nullX = ReferenceEquals(null, x), nullY = ReferenceEquals(null, y);
return nullX || nullY ? nullX && nullY : x.Equals(y);
}

public static bool operator !=(Onion x, Onion y) { return !(x == y); }

public static Onion operator +(Onion x, Onion y) { return Operate(x, y, (a, b) => a + b); }

public static Onion operator -(Onion x, Onion y) { return Operate(x, y, (a, b) => a - b); }

public static Onion operator *(Onion x, Onion y)
{
int dimX = x.Dim, dimY = y.Dim;
if (dimX == 1) return dimY == 1 ? x[0]*y[0] : Operate(y, value => x[0]*value);
if (dimY == 1) return Operate(x, value => value*y[0]);
if (dimX < dimY) return new Onion(x*y.A, y.B*x);
if (dimX > dimY) return new Onion(x.A*y, x.B*~y);
Onion a = x.A, b = x.B, c = y.A, d = y.B;
return new Onion(a*c - ~d*b, d*a + b*~c);
}

public static Onion operator /(Onion x, Onion y) { return x*y.Reciprocal; }

public static Onion operator ^(Onion x, long y)
{
if (y < 0) return x.Reciprocal ^ -y;
Onion z = 1;
var power = x;
var unity = true;
for (long mask = 1; y > 0 && mask > 0; mask <<= 1)
{
if ((y & mask) != 0)
{
if (unity) { z = power; unity = false; }
else z *= power;
y &= ~mask;
}
if (mask > 0) power *= power;
}
return z;
}

public static Onion operator ^(Onion x, Onion y)
{
return Exp(Log(x)*y);
}

private static void AppendDim(StringBuilder result, int index)
{
for (int n = 0, keys = index, mask = 1; keys != 0 && mask > 0; n++, mask <<= 1)
if ((keys & mask) != 0)
{
result.Append((char)('i' + n));
keys &= ~mask;
}
}

public override bool Equals(object obj)
{
return
!ReferenceEquals(null, obj) &&
(ReferenceEquals(this, obj) || obj is Onion && Equals((Onion) obj));
}

private bool Equals(Onion other)
{
if (ReferenceEquals(null, other)) return false;
if (ReferenceEquals(this, other)) return true;
for (var n = 0; n < Math.Max(Dim, other.Dim); n++)
if (this[n] != other[n]) return false;
return true;
}

public static Onion Exp(Onion x)
{
var re = x.Re;
Onion z = Math.Exp(re);
if (x.Dim > 1)
{
var im = x.Im;
var angle = im.Norm;
if (angle != 0) z *= (Math.Cos(angle) + im*(Math.Sin(angle)/angle));
}
return z;
}

public override int GetHashCode()
{
return Values
.Select(value => value.GetHashCode())
.Aggregate((u, v) => u ^ v);
}

public static Onion Log(Onion x)
{
var norm = x.Norm;
Onion z = Math.Log(norm);
if (x.Dim > 1)
{
var im = x.Im;
var angle = im.Norm;
if (angle != 0) z += im*(Math.Acos(x.Re/norm)/angle);
}
return z;
}

private static Onion Operate(Onion x, Func<double, double> f)
{
var z = new Onion(x);
for (var n = 0; n < x.Dim; n++) z.Values[n] = f(z[n]);
return z;
}

private static Onion Operate(Onion x, Onion y, Func<double, double, double> f)
{
var dim = Math.Max(x.Dim, y.Dim);
var z = new Onion { Dim = dim };
for (var n = 0; n < dim; n++) z.Values[n] = f(x[n], y[n]);
return z;
}

public override string ToString()
{
var result = new StringBuilder();
var empty = true;
for (var n = 0; n < Dim; n++)
{
var value = this[n];
if (value == 0) continue;
if (empty)
{
if (value < 0) result.Append('-');
empty = false;
}
else result.Append(value < 0 ? " - " : " + ");
var absValue = Math.Abs(value);
if (n == 0 || absValue != 1) result.Append(absValue);
AppendDim(result, n);
}
if (empty) result.Append('0');
return result.ToString();
}
}
}
Happy hypercomplex numbering.

No comments:

Post a Comment