From: Michael Orlitzky <michael@orlitzky.com>
Date: Thu, 25 Feb 2021 01:31:48 +0000 (-0500)
Subject: eja: make gram_schmidt work in Cartesian product algebras.
X-Git-Url: http://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=8492b6f01186a5c5df9f7ccb007604958ffdf1f1;hp=a0e9fb9e4c56781b9dad824bac517ca3e56ef461;p=sage.d.git

eja: make gram_schmidt work in Cartesian product algebras.
---

diff --git a/mjo/eja/eja_utils.py b/mjo/eja/eja_utils.py
index 29edf5b..e8ed4db 100644
--- a/mjo/eja/eja_utils.py
+++ b/mjo/eja/eja_utils.py
@@ -110,9 +110,6 @@ def gram_schmidt(v, inner_product=None):
         inner_product = lambda x,y: x.inner_product(y)
     norm = lambda x: inner_product(x,x).sqrt()
 
-    def proj(x,y):
-        return (inner_product(x,y)/inner_product(x,x))*x
-
     v = list(v) # make a copy, don't clobber the input
 
     # Drop all zero vectors before we start.
@@ -124,10 +121,26 @@ def gram_schmidt(v, inner_product=None):
 
     R = v[0].base_ring()
 
+    # Define a scaling operation that can be used on tuples.
+    # Oh and our "zero" needs to belong to the right space.
+    scale = lambda x,alpha: x*alpha
+    zero = v[0].parent().zero()
+    if hasattr(v[0], 'cartesian_factors'):
+        P = v[0].parent()
+        scale = lambda x,alpha: P(tuple( x_i*alpha
+                                         for x_i in x.cartesian_factors() ))
+
+
+    def proj(x,y):
+        return scale(x, (inner_product(x,y)/inner_product(x,x)))
+
     # First orthogonalize...
     for i in range(1,len(v)):
         # Earlier vectors can be made into zero so we have to ignore them.
-        v[i] -= sum( proj(v[j],v[i]) for j in range(i) if not v[j].is_zero() )
+        v[i] -= sum( (proj(v[j],v[i])
+                      for j in range(i)
+                      if not v[j].is_zero() ),
+                     zero )
 
     # And now drop all zero vectors again if they were "orthogonalized out."
     v = [ v_i for v_i in v if not v_i.is_zero() ]
@@ -136,6 +149,6 @@ def gram_schmidt(v, inner_product=None):
     # them here because then our subalgebra would have a bigger field
     # than the superalgebra.
     for i in range(len(v)):
-        v[i] = v[i] / norm(v[i])
+        v[i] = scale(v[i], ~norm(v[i]))
 
     return v