An earlier post included code for multiplying quaternions, octonions, and sedenions. The code was a little clunky, so I refactor it here.
def conj(x): xstar = -x xstar *= -1 return xstar def CayleyDickson(x, y): n = len(x) if n == 1: return x*y m = n // 2 a, b = x[:m], x[m:] c, d = y[:m], y[m:] z = np.zeros(n) z[:m] = CayleyDickson(a, c) - CayleyDickson(conj(d), b) z[m:] = CayleyDickson(d, a) + CayleyDickson(b, conj(c)) return z
CayleyDickson function implements the Cayley-Dickson construction and can be used to multiply real, complex, quaternion, and octonion numbers. In fact, it can be used to implement multiplication in any real vector space of dimension 2n. The numeric types listed above correspond to n = 0, 1, 2, and 3. These are the only normed division algebras over the reals.
When n = 4 we get the sedenions, which are not a division algebra because they contain zero divisors, and the code can be used for any larger value of n as well. As noted before, the algebraic properties degrade as n increases, though I don’t think they get any worse after n = 4.
If you wanted to make the code more robust, you could add code to verify that the arguments
y have the same length, and that their common length is a power of 2. (You could just check that the length is either 1 or even; if it’s not a power of 2 the recursion will eventually produce an odd argument.)