Overview
--------
This is a collection of design notes that should eventually wind up in
the documentation. I'm just not sure where they go yet.

Matrix representations
----------------------

All algebras have a "matrix representation" of their elements. This is
the original, ambient representation of the elements as either column
(n-by-1) or square (n-by-n) matrices. For example, the elements of the
Jordan spin algebra are column vectors, and the elements of real
symmetric matrix algebra are... real symmetric matrices.

The CombinatorialFreeModule class actually supports such an
alternative representation of its generators since it subclasses
IndexedGenerators. However, using matrices as the index set turns out
to be ugly: printing the generators, and especially printing an
algebra element as a sum of said generators comes out wonky, since
the matrices require more than one line. For example,

    sage: A = matrix(QQ,[[1,2],[3,4]])
    sage: B = matrix(QQ,[[5,6],[7,8]])
    sage: A.set_immutable()
    sage: B.set_immutable()
    sage: M = CombinatorialFreeModule(QQ,[A,B],bracket=False,prefix="")
    sage: 2*M(A) + 3*M(B)
    2*[1 2]
    [3 4] + 3*[5 6]
    [7 8]

And things only get worse if you leave the prefix in there to
distinguish between e.g. the super- and sub-algebra elements
corresponding to the same matrix. Thus, we store out own copy
of the matrix generators, and have our own set of methods for
accessing them.

Why allow matrix representations for all algebras, rather than just
the matrix algebras?

1. We already have a to_vector() operation that turns an algebra
   element into a vector whose coordinates live in the algebra's
   base_ring(). Adding a to_matrix() operation is a natural
   generalization of that.

2. Having access to the ambient coordinates in a general way is useful
   when converting between other coordinate systems. If we have two
   subalgebras B and C of A, we can use to_matrix() to go from, say,
   C -> A -> B rather than having to convert from C to B directly.


Fielding questions
------------------

All Euclidean Jordan algebras are over the real scalar field. This
presents a problem: in SageMath, the matrix and vector classes don't
support scalar fields that are different than their entries. And three
of the simple families of Euclidean Jordan algebras are matrices with
non-real entries: the Hermitian comples, quaternion, and octonion
algebras.

At least in the complex and quaternion case, we can "embed" the
complex numbers and quaternions into the space of 2-by-2 or 4-by-4
matrices. But the octonions are not associative, so they can't be
embedded (via a homomorphism) into any real matrix space. So what
do we do? Write it ourselves, obviously.

In contrast to the algebra of real symmetric matrices, the complex,
quaternion, and octonion matrix algebras are implemented separately,
as a subclasses of CombinatorialFreeModule, to work around that
issue. The custom class supports a scalar field that is different than
the entries of the matrices. However, this means that we actually have
FOUR different types of "matrices" to support:

  (1) Sage vectors
  (2) Sage matrices
  (3) Our custom matrices
  (4) Cartesian products of the (1) through (3)

The real symmetric matrices could of course be implemented in the same
manner as the others; but for the sake of the user interface, we must
also support at least the usual SageMath vectors and matrices. Having
the real symmetric matrices actually be (SageMath) matrices ensures
that we don't accidentally break support for such things.

Note: this has one less-than-obvious consequence: we have to assume
that the user has supplied an entirely-correct basis (with entries in
the correct structure). We generally cannot mess witht the entries of
his basis, or use them to figure out what (for example) the ambient
scalar ring is. None of these are insurmountable obstacles; we just
have to be a little careful distinguishing between what's inside the
algebra elements and what's outside them.

Basis normalization
-------------------
For performance reasons, we prefer the algebra constructor to
orthonormalize its own basis. We _could_ ask the user to do that,
but there's a good reason to do it ourselves: if _we_ run
Gram-Schmidt, then we can compute/store the matrix that undoes the
process. Undoing the change-of-coordinates allows us to perform some
computations in the original basis (like the "characteristic
polynomial of"), which can be much faster when the original basis
contains only rational entries.

Debugging
---------
There are several places in the code where we set check_field=False
and check_axioms=False because the theory guarantees that they will be
satisfied. Well, you know what they say about theory and practice. The
first thing you should do when a problem is discovered it replace all
of those with check_field=True and check_axioms=True, and then re-run
the test suite. The Cartesian product class bypasses its superclass
constructor, so any extra axiom/field checks on product algebras must
be inserted at debug-time.