--- /dev/null
+from sage.misc.cachefunc import cached_method
+from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
+from sage.combinat.free_module import CombinatorialFreeModule
+from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement
+from sage.categories.magmatic_algebras import MagmaticAlgebras
+from sage.rings.all import AA, ZZ
+from sage.matrix.matrix_space import MatrixSpace
+from sage.misc.table import table
+
+from mjo.matrix_algebra import MatrixAlgebra, MatrixAlgebraElement
+
+class Octonion(IndexedFreeModuleElement):
+ def conjugate(self):
+ r"""
+ SETUP::
+
+ sage: from mjo.hurwitz import Octonions
+
+ EXAMPLES::
+
+ sage: O = Octonions()
+ sage: x = sum(O.gens())
+ sage: x.conjugate()
+ e0 - e1 - e2 - e3 - e4 - e5 - e6 - e7
+
+ TESTS::
+
+ Conjugating twice gets you the original element::
+
+ sage: set_random_seed()
+ sage: O = Octonions()
+ sage: x = O.random_element()
+ sage: x.conjugate().conjugate() == x
+ True
+
+ """
+ C = MatrixSpace(ZZ,8).diagonal_matrix((1,-1,-1,-1,-1,-1,-1,-1))
+ return self.parent().from_vector(C*self.to_vector())
+
+ def real(self):
+ r"""
+ Return the real part of this octonion.
+
+ The real part of an octonion is its projection onto the span
+ of the first generator. In other words, the "first dimension"
+ is real and the others are imaginary.
+
+ SETUP::
+
+ sage: from mjo.hurwitz import Octonions
+
+ EXAMPLES::
+
+ sage: O = Octonions()
+ sage: x = sum(O.gens())
+ sage: x.real()
+ e0
+
+ TESTS:
+
+ This method is idempotent::
+
+ sage: set_random_seed()
+ sage: O = Octonions()
+ sage: x = O.random_element()
+ sage: x.real().real() == x.real()
+ True
+
+ """
+ return (self + self.conjugate())/2
+
+ def imag(self):
+ r"""
+ Return the imaginary part of this octonion.
+
+ The imaginary part of an octonion is its projection onto the
+ orthogonal complement of the span of the first generator. In
+ other words, the "first dimension" is real and the others are
+ imaginary.
+
+ SETUP::
+
+ sage: from mjo.hurwitz import Octonions
+
+ EXAMPLES::
+
+ sage: O = Octonions()
+ sage: x = sum(O.gens())
+ sage: x.imag()
+ e1 + e2 + e3 + e4 + e5 + e6 + e7
+
+ TESTS:
+
+ This method is idempotent::
+
+ sage: set_random_seed()
+ sage: O = Octonions()
+ sage: x = O.random_element()
+ sage: x.imag().imag() == x.imag()
+ True
+
+ """
+ return (self - self.conjugate())/2
+
+ def _norm_squared(self):
+ return (self*self.conjugate()).coefficient(0)
+
+ def norm(self):
+ r"""
+ Return the norm of this octonion.
+
+ SETUP::
+
+ sage: from mjo.hurwitz import Octonions
+
+ EXAMPLES::
+
+ sage: O = Octonions()
+ sage: O.one().norm()
+ 1
+
+ TESTS:
+
+ The norm is nonnegative and belongs to the base field::
+
+ sage: set_random_seed()
+ sage: O = Octonions()
+ sage: n = O.random_element().norm()
+ sage: n >= 0 and n in O.base_ring()
+ True
+
+ The norm is homogeneous::
+
+ sage: set_random_seed()
+ sage: O = Octonions()
+ sage: x = O.random_element()
+ sage: alpha = O.base_ring().random_element()
+ sage: (alpha*x).norm() == alpha.abs()*x.norm()
+ True
+
+ """
+ return self._norm_squared().sqrt()
+
+ # The absolute value notation is typically used for complex numbers...
+ # and norm() isn't supported in AA, so this lets us use abs() in all
+ # of the division algebras we need.
+ abs = norm
+
+ def inverse(self):
+ r"""
+ Return the inverse of this element if it exists.
+
+ SETUP::
+
+ sage: from mjo.hurwitz import Octonions
+
+ EXAMPLES::
+
+ sage: O = Octonions()
+ sage: x = sum(O.gens())
+ sage: x*x.inverse() == O.one()
+ True
+
+ ::
+
+ sage: O = Octonions()
+ sage: O.one().inverse() == O.one()
+ True
+
+ TESTS::
+
+ sage: set_random_seed()
+ sage: O = Octonions()
+ sage: x = O.random_element()
+ sage: x.is_zero() or ( x*x.inverse() == O.one() )
+ True
+
+ """
+ if self.is_zero():
+ raise ValueError("zero is not invertible")
+ return self.conjugate()/self._norm_squared()
+
+
+ def cayley_dickson(self, Q=None):
+ r"""
+ Return the Cayley-Dickson representation of this element in terms
+ of the quaternion algebra ``Q``.
+
+ The Cayley-Dickson representation is an identification of
+ octionions `x` and `y` with pairs of quaternions `(a,b)` and
+ `(c,d)` respectively such that:
+
+ * `x + y = (a+b, c+d)`
+ * `xy` = (ac - \bar{d}*b, da + b\bar{c})`
+ * `\bar{x} = (a,-b)`
+
+ where `\bar{x}` denotes the conjugate of `x`.
+
+ SETUP::
+
+ sage: from mjo.hurwitz import Octonions
+
+ EXAMPLES::
+
+ sage: O = Octonions()
+ sage: x = sum(O.gens())
+ sage: x.cayley_dickson()
+ (1 + i + j + k, 1 + i + j + k)
+
+ """
+ if Q is None:
+ Q = QuaternionAlgebra(self.base_ring(), -1, -1)
+
+ i,j,k = Q.gens()
+ a = (self.coefficient(0)*Q.one() +
+ self.coefficient(1)*i +
+ self.coefficient(2)*j +
+ self.coefficient(3)*k )
+ b = (self.coefficient(4)*Q.one() +
+ self.coefficient(5)*i +
+ self.coefficient(6)*j +
+ self.coefficient(7)*k )
+
+ from sage.categories.sets_cat import cartesian_product
+ P = cartesian_product([Q,Q])
+ return P((a,b))
+
+
+class Octonions(CombinatorialFreeModule):
+ r"""
+ SETUP::
+
+ sage: from mjo.hurwitz import Octonions
+
+ EXAMPLES::
+
+ sage: Octonions()
+ Octonion algebra with base ring Algebraic Real Field
+ sage: Octonions(field=QQ)
+ Octonion algebra with base ring Rational Field
+
+ """
+ def __init__(self,
+ field=AA,
+ prefix="e"):
+
+ # Not associative, not commutative
+ category = MagmaticAlgebras(field).FiniteDimensional()
+ category = category.WithBasis().Unital()
+
+ super().__init__(field,
+ range(8),
+ element_class=Octonion,
+ category=category,
+ prefix=prefix,
+ bracket=False)
+
+ # The product of each basis element is plus/minus another
+ # basis element that can simply be looked up on
+ # https://en.wikipedia.org/wiki/Octonion
+ e0, e1, e2, e3, e4, e5, e6, e7 = self.gens()
+ self._multiplication_table = (
+ (e0, e1, e2, e3, e4, e5, e6, e7),
+ (e1,-e0, e3,-e2, e5,-e4,-e7, e6),
+ (e2,-e3,-e0, e1, e6, e7,-e4,-e5),
+ (e3, e2,-e1,-e0, e7,-e6, e5,-e4),
+ (e4,-e5,-e6,-e7,-e0, e1, e2, e3),
+ (e5, e4,-e7, e6,-e1,-e0,-e3, e2),
+ (e6, e7, e4,-e5,-e2, e3,-e0,-e1),
+ (e7,-e6, e5, e4,-e3,-e2, e1,-e0),
+ )
+
+ def product_on_basis(self, i, j):
+ return self._multiplication_table[i][j]
+
+ def one_basis(self):
+ r"""
+ Return the monomial index (basis element) corresponding to the
+ octonion unit element.
+
+ SETUP::
+
+ sage: from mjo.hurwitz import Octonions
+
+ TESTS:
+
+ This gives the correct unit element::
+
+ sage: set_random_seed()
+ sage: O = Octonions()
+ sage: x = O.random_element()
+ sage: x*O.one() == x and O.one()*x == x
+ True
+
+ """
+ return 0
+
+ def _repr_(self):
+ return ("Octonion algebra with base ring %s" % self.base_ring())
+
+ def multiplication_table(self):
+ """
+ Return a visual representation of this algebra's multiplication
+ table (on basis elements).
+
+ SETUP::
+
+ sage: from mjo.hurwitz import Octonions
+
+ EXAMPLES:
+
+ The multiplication table is what Wikipedia says it is::
+
+ sage: Octonions().multiplication_table()
+ +----++----+-----+-----+-----+-----+-----+-----+-----+
+ | * || e0 | e1 | e2 | e3 | e4 | e5 | e6 | e7 |
+ +====++====+=====+=====+=====+=====+=====+=====+=====+
+ | e0 || e0 | e1 | e2 | e3 | e4 | e5 | e6 | e7 |
+ +----++----+-----+-----+-----+-----+-----+-----+-----+
+ | e1 || e1 | -e0 | e3 | -e2 | e5 | -e4 | -e7 | e6 |
+ +----++----+-----+-----+-----+-----+-----+-----+-----+
+ | e2 || e2 | -e3 | -e0 | e1 | e6 | e7 | -e4 | -e5 |
+ +----++----+-----+-----+-----+-----+-----+-----+-----+
+ | e3 || e3 | e2 | -e1 | -e0 | e7 | -e6 | e5 | -e4 |
+ +----++----+-----+-----+-----+-----+-----+-----+-----+
+ | e4 || e4 | -e5 | -e6 | -e7 | -e0 | e1 | e2 | e3 |
+ +----++----+-----+-----+-----+-----+-----+-----+-----+
+ | e5 || e5 | e4 | -e7 | e6 | -e1 | -e0 | -e3 | e2 |
+ +----++----+-----+-----+-----+-----+-----+-----+-----+
+ | e6 || e6 | e7 | e4 | -e5 | -e2 | e3 | -e0 | -e1 |
+ +----++----+-----+-----+-----+-----+-----+-----+-----+
+ | e7 || e7 | -e6 | e5 | e4 | -e3 | -e2 | e1 | -e0 |
+ +----++----+-----+-----+-----+-----+-----+-----+-----+
+
+ """
+ n = self.dimension()
+ # Prepend the header row.
+ M = [["*"] + list(self.gens())]
+
+ # And to each subsequent row, prepend an entry that belongs to
+ # the left-side "header column."
+ M += [ [self.monomial(i)] + [ self.monomial(i)*self.monomial(j)
+ for j in range(n) ]
+ for i in range(n) ]
+
+ return table(M, header_row=True, header_column=True, frame=True)
+
+
+
+
+
+class HurwitzMatrixAlgebraElement(MatrixAlgebraElement):
+ def is_hermitian(self):
+ r"""
+
+ SETUP::
+
+ sage: from mjo.hurwitz import HurwitzMatrixAlgebra
+
+ EXAMPLES::
+
+ sage: A = HurwitzMatrixAlgebra(QQbar, ZZ, 2)
+ sage: M = A([ [ 0,I],
+ ....: [-I,0] ])
+ sage: M.is_hermitian()
+ True
+
+ """
+ return all( self[i,j] == self[j,i].conjugate()
+ for i in range(self.nrows())
+ for j in range(self.ncols()) )
+
+
+class HurwitzMatrixAlgebra(MatrixAlgebra):
+ r"""
+ A class of matrix algebras whose entries come from a Hurwitz
+ algebra.
+
+ For our purposes, we consider "a Hurwitz" algebra to be the real
+ or complex numbers, the quaternions, or the octonions. These are
+ typically also referred to as the Euclidean Hurwitz algebras, or
+ the normed division algebras.
+
+ By the Cayley-Dickson construction, each Hurwitz algebra is an
+ algebra over the real numbers, so we restrict the scalar field in
+ this case to be real. This also allows us to more accurately
+ produce the generators of the matrix algebra.
+ """
+ Element = HurwitzMatrixAlgebraElement
+
+ def __init__(self, entry_algebra, scalars, n, **kwargs):
+ from sage.rings.all import RR
+ if not scalars.is_subring(RR):
+ # Not perfect, but it's what we're using.
+ raise ValueError("scalar field is not real")
+
+ super().__init__(entry_algebra, scalars, n, **kwargs)
+
+ def entry_algebra_gens(self):
+ r"""
+ Return the generators of (that is, a basis for) the entries of
+ this matrix algebra.
+
+ This works around the inconsistency in the ``gens()`` methods
+ of the real/complex numbers, quaternions, and octonions.
+
+ SETUP::
+
+ sage: from mjo.hurwitz import Octonions, HurwitzMatrixAlgebra
+
+ EXAMPLES:
+
+ The inclusion of the unit element is inconsistent across
+ (subalgebras of) Hurwitz algebras::
+
+ sage: AA.gens()
+ (1,)
+ sage: QQbar.gens()
+ (I,)
+ sage: QuaternionAlgebra(AA,1,-1).gens()
+ [i, j, k]
+ sage: Octonions().gens()
+ (e0, e1, e2, e3, e4, e5, e6, e7)
+
+ The unit element is always returned by this method, so the
+ sets of generators have cartinality 1,2,4, and 8 as you'd
+ expect::
+
+ sage: HurwitzMatrixAlgebra(AA, AA, 2).entry_algebra_gens()
+ (1,)
+ sage: HurwitzMatrixAlgebra(QQbar, AA, 2).entry_algebra_gens()
+ (1, I)
+ sage: Q = QuaternionAlgebra(AA,-1,-1)
+ sage: HurwitzMatrixAlgebra(Q, AA, 2).entry_algebra_gens()
+ (1, i, j, k)
+ sage: O = Octonions()
+ sage: HurwitzMatrixAlgebra(O, AA, 2).entry_algebra_gens()
+ (e0, e1, e2, e3, e4, e5, e6, e7)
+
+ """
+ gs = self.entry_algebra().gens()
+ one = self.entry_algebra().one()
+ if one in gs:
+ return gs
+ else:
+ return (one,) + tuple(gs)
+
+
+
+class OctonionMatrixAlgebra(HurwitzMatrixAlgebra):
+ r"""
+ The algebra of ``n``-by-``n`` matrices with octonion entries over
+ (a subfield of) the real numbers.
+
+ The usual matrix spaces in SageMath don't support octonion entries
+ because they assume that the entries of the matrix come from a
+ commutative and associative ring, and the octonions are neither.
+
+ SETUP::
+
+ sage: from mjo.hurwitz import OctonionMatrixAlgebra
+
+ EXAMPLES::
+
+ sage: OctonionMatrixAlgebra(3)
+ Module of 3 by 3 matrices with entries in Octonion algebra with base
+ ring Algebraic Real Field over the scalar ring Algebraic Real Field
+ sage: OctonionMatrixAlgebra(3,QQ)
+ Module of 3 by 3 matrices with entries in Octonion algebra with base
+ ring Rational Field over the scalar ring Rational Field
+
+ ::
+
+ sage: A = OctonionMatrixAlgebra(2)
+ sage: e0,e1,e2,e3,e4,e5,e6,e7 = A.entry_algebra().gens()
+ sage: A([ [e0+e4, e1+e5],
+ ....: [e2-e6, e3-e7] ])
+ +---------+---------+
+ | e0 + e4 | e1 + e5 |
+ +---------+---------+
+ | e2 - e6 | e3 - e7 |
+ +---------+---------+
+
+ ::
+
+ sage: A1 = OctonionMatrixAlgebra(1,QQ)
+ sage: A2 = OctonionMatrixAlgebra(1,QQ)
+ sage: cartesian_product([A1,A2])
+ Module of 1 by 1 matrices with entries in Octonion algebra with
+ base ring Rational Field over the scalar ring Rational Field (+)
+ Module of 1 by 1 matrices with entries in Octonion algebra with
+ base ring Rational Field over the scalar ring Rational Field
+
+ TESTS::
+
+ sage: set_random_seed()
+ sage: A = OctonionMatrixAlgebra(ZZ.random_element(10))
+ sage: x = A.random_element()
+ sage: x*A.one() == x and A.one()*x == x
+ True
+
+ """
+ def __init__(self, n, scalars=AA, prefix="E", **kwargs):
+ super().__init__(Octonions(field=scalars),
+ scalars,
+ n,
+ prefix=prefix,
+ **kwargs)
+
+class QuaternionMatrixAlgebra(HurwitzMatrixAlgebra):
+ r"""
+ The algebra of ``n``-by-``n`` matrices with quaternion entries over
+ (a subfield of) the real numbers.
+
+ The usual matrix spaces in SageMath don't support quaternion entries
+ because they assume that the entries of the matrix come from a
+ commutative ring, and the quaternions are not commutative.
+
+ SETUP::
+
+ sage: from mjo.hurwitz import QuaternionMatrixAlgebra
+
+ EXAMPLES::
+
+ sage: QuaternionMatrixAlgebra(3)
+ Module of 3 by 3 matrices with entries in Quaternion
+ Algebra (-1, -1) with base ring Algebraic Real Field
+ over the scalar ring Algebraic Real Field
+ sage: QuaternionMatrixAlgebra(3,QQ)
+ Module of 3 by 3 matrices with entries in Quaternion
+ Algebra (-1, -1) with base ring Rational Field over
+ the scalar ring Rational Field
+
+ ::
+
+ sage: A = QuaternionMatrixAlgebra(2)
+ sage: i,j,k = A.entry_algebra().gens()
+ sage: A([ [1+i, j-2],
+ ....: [k, k+j] ])
+ +-------+--------+
+ | 1 + i | -2 + j |
+ +-------+--------+
+ | k | j + k |
+ +-------+--------+
+
+ ::
+
+ sage: A1 = QuaternionMatrixAlgebra(1,QQ)
+ sage: A2 = QuaternionMatrixAlgebra(2,QQ)
+ sage: cartesian_product([A1,A2])
+ Module of 1 by 1 matrices with entries in Quaternion Algebra
+ (-1, -1) with base ring Rational Field over the scalar ring
+ Rational Field (+) Module of 2 by 2 matrices with entries in
+ Quaternion Algebra (-1, -1) with base ring Rational Field over
+ the scalar ring Rational Field
+
+ TESTS::
+
+ sage: set_random_seed()
+ sage: A = QuaternionMatrixAlgebra(ZZ.random_element(10))
+ sage: x = A.random_element()
+ sage: x*A.one() == x and A.one()*x == x
+ True
+
+ """
+ def __init__(self, n, scalars=AA, **kwargs):
+ # The -1,-1 gives us the "usual" definition of quaternion
+ Q = QuaternionAlgebra(scalars,-1,-1)
+ super().__init__(Q, scalars, n, **kwargs)