From c9db5fa5d2f30833e4da85d32418048d34483e43 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Fri, 27 Nov 2020 10:38:18 -0500
Subject: [PATCH] eja: use new constructor for BilinearFormEJA.

---
 mjo/eja/eja_algebra.py | 75 ++++++++++++++++++++++++------------------
 1 file changed, 43 insertions(+), 32 deletions(-)

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index 4609ca2..9eba669 100644
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -1062,8 +1062,15 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge
 
         if orthonormalize:
             from mjo.eja.eja_utils import gram_schmidt
+            A = matrix(field, vector_basis)
             vector_basis = gram_schmidt(vector_basis, inner_product)
             W = V.span_of_basis( vector_basis )
+            Q = matrix(field, vector_basis)
+            # A = QR <==> A.T == R.T*Q.T
+            # So, Q.solve_right() is equivalent to the Q.T.solve_left()
+            # that we want.
+            self._deorthonormalization_matrix = Q.solve_right(A)
+
             if basis_is_matrices:
                 from mjo.eja.eja_utils import _vec2mat
                 basis = tuple( map(_vec2mat,vector_basis) )
@@ -1073,15 +1080,30 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge
         mult_table = [ [0 for i in range(n)] for j in range(n) ]
         ip_table = [ [0 for i in range(n)] for j in range(n) ]
 
+        # Note: the Jordan and inner- products are defined in terms
+        # of the ambient basis. It's important that their arguments
+        # are in ambient coordinates as well.
         for i in range(n):
             for j in range(i+1):
-                # do another mat2vec because the multiplication
-                # table is in terms of vectors
-                elt = _mat2vec(jordan_product(basis[i],basis[j]))
+                # ortho basis w.r.t. ambient coords
+                q_i = vector_basis[i]
+                q_j = vector_basis[j]
+
+                if basis_is_matrices:
+                    q_i = _vec2mat(q_i)
+                    q_j = _vec2mat(q_j)
+
+                elt = jordan_product(q_i, q_j)
+                ip = inner_product(q_i, q_j)
+
+                if basis_is_matrices:
+                    # do another mat2vec because the multiplication
+                    # table is in terms of vectors
+                    elt = _mat2vec(elt)
+
                 elt = W.coordinate_vector(elt)
                 mult_table[i][j] = elt
                 mult_table[j][i] = elt
-                ip = inner_product(basis[i],basis[j])
                 ip_table[i][j] = ip
                 ip_table[j][i] = ip
 
@@ -2193,7 +2215,7 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebraNg,
         return cls(n, field, **kwargs)
 
 
-class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
+class BilinearFormEJA(RationalBasisEuclideanJordanAlgebraNg,
                       ConcreteEuclideanJordanAlgebra):
     r"""
     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
@@ -2274,39 +2296,28 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
         True
     """
     def __init__(self, B, field=AA, **kwargs):
-        n = B.nrows()
-
         if not B.is_positive_definite():
             raise ValueError("bilinear form is not positive-definite")
 
+        n = B.nrows()
         V = VectorSpace(field, n)
-        mult_table = [[V.zero() for j in range(n)] for i in range(n)]
-        for i in range(n):
-            for j in range(n):
-                x = V.gen(i)
-                y = V.gen(j)
-                x0 = x[0]
-                xbar = x[1:]
-                y0 = y[0]
-                ybar = y[1:]
-                z0 = (B*x).inner_product(y)
-                zbar = y0*xbar + x0*ybar
-                z = V([z0] + zbar.list())
-                mult_table[i][j] = z
-
-        # Inner products are real numbers and not algebra
-        # elements, so once we turn the algebra element
-        # into a vector in inner_product(), we never go
-        # back. As a result -- contrary to what we do with
-        # self._multiplication_table -- we store the inner
-        # product table as a plain old matrix and not as
-        # an algebra operator.
-        ip_table = B
-        self._inner_product_matrix = ip_table
+
+        def inner_product(x,y):
+            return (B*x).inner_product(y)
+
+        def jordan_product(x,y):
+            x0 = x[0]
+            xbar = x[1:]
+            y0 = y[0]
+            ybar = y[1:]
+            z0 = inner_product(x,y)
+            zbar = y0*xbar + x0*ybar
+            return V([z0] + zbar.list())
 
         super(BilinearFormEJA, self).__init__(field,
-                                              mult_table,
-                                              check_axioms=False,
+                                              V.basis(),
+                                              jordan_product,
+                                              inner_product,
                                               **kwargs)
 
         # The rank of this algebra is two, unless we're in a
-- 
2.44.2