# A,B - algebras with linear maps f:A->A and g:B->B; 
# returns the pair (A\otimes B, the linear map f\otimes g)
TensorProductOfAlgebrasWithMaps := function (A, f, B, g)
    local AotimesB, basisA, basisB, basisAotimesB, ca, cb, dimA, dimB, 
          i, images, j, K, k, l, list, m, matrix, n, row, scA, scB, T;

    K := LeftActingDomain(A);
    if (K <> LeftActingDomain(B)) then
        Error ("tensor factors should be defined over the same field");
    fi;

    dimA := Dimension (A);
    dimB := Dimension (B);
    basisA := CanonicalBasis (A);
    basisB := CanonicalBasis (B);
    scA := StructureConstantsTable (basisA);
    scB := StructureConstantsTable (basisB);

    # structure constants of the resulting algebra
    T := EmptySCTable (dimA * dimB, Zero(K));  

    # list of images with respect to T\otimes S of basic elements in A\otimes B
    images := [];

    # matrix of the linear map f \otimes g
    matrix := [];

    # assuming (a_i) is a basis of A and (b_i) is a basis B,
    # basis of A\otimes B is a set of pairs (a_i,b_j) numbered in the reverse
    # lexicographic order (so (a_i,b_j) is at i + (j-1)*Dimension(A) place)
    for j in [1 .. dimB] do
        cb := Coefficients (basisB, Image (g, basisB[j]));
        for i in [1 .. dimA] do
            ca := Coefficients (basisA, Image (f, basisA[i]));

            # the row to be added to the matrix
            row := [];

            for l in [1 .. dimB] do
                for k in [1 .. dimA] do
                    list := [];  # auxiliary list used to form structure 
                                 # constants of the resulting algebra
                    for n in [1 .. Length (scA[i][k][1])] do
                        for m in [1 .. Length (scB[j][l][1])] do
                            Add (list, scA[i][k][2][n] * scB[j][l][2][m]);
                            Add (list, 
                                scA[i][k][1][n] + (scB[j][l][1][m] - 1)*dimA);
                        od;
                    od;
                    SetEntrySCTable (T, i+(j-1)*dimA, k+(l-1)*dimA, list);
                    Add (row, ca[k]*cb[l]);
                od;
            od;
            Add (matrix, row);
        od;
    od;
    
    AotimesB := AlgebraByStructureConstants (K, T);
    basisAotimesB := CanonicalBasis (AotimesB);

    return ([AotimesB, LeftModuleHomomorphismByMatrix 
                (basisAotimesB, matrix, basisAotimesB)]);
end;


# given two vector spaces V, W of the same dimension, and with canononical bases
# and a linear map f:V->V, produces the linear map W->W with the same matrix
# as f
Map2Map := function (V, W, f)
    local b, bV, bW, images;
    bV := CanonicalBasis (V);
    bW := CanonicalBasis (W);
    images := [];
    for b in bV do
        Add (images, LinearCombination (bW, Coefficients (bV, Image (f, b))));
    od;
    return (LeftModuleHomomorphismByImages (W, W, bW, images));
end;



# transposition on the matrix algebra A = M_n(K)
Transposition := function (A)
    local b;
    b := CanonicalBasis (A);
    return (LeftModuleHomomorphismByImages (
        A, A, b, List (b, b -> TransposedMat (b))));
end;



# for an algebra A with involution J, constructs the algebra S^-(A,J) of
# skew-symmetric elements with respect to the involution
SminusAJ := function (A, J)
    local Aminus, Jminus, gen, b;
    Aminus := MinusAlgebra (A);
    Jminus := Map2Map (A, Aminus, J);
    gen := [];
    for b in CanonicalBasis (Aminus) do
        Add (gen, b - Image (Jminus, b));
    od;
    return (Subalgebra (Aminus, gen));
end;



# for an algebra A with involution J, constructs the algebra S^+(A,J) of
# symmetric elements with respect to the involution
SplusAJ := function (A, J)
    local Aplus, Jplus, gen, b;
    Aplus := PlusAlgebra (A);
    Jplus := Map2Map (A, Aplus, J);
    gen := [];
    for b in CanonicalBasis (Aplus) do
        Add (gen, b + Image (Jplus, b));
    od;
    return (Subalgebra (Aplus, gen));
end;



# returns M_n(O) with its standard involution
# K: ground field
MatrixOctonionAlgebra := function (K, n)
    local M, O;
    M := MatrixAlgebra (K, n);
    O := OctonionAlgebra (K, -1, -1, -1);
    return (TensorProductOfAlgebrasWithMaps 
                (M, Transposition (M), O, OctonionInvolution (O)));
end;



HermitianMatrixOctonionAlgebra := function (K, n)
    local m;
    m := MatrixOctonionAlgebra (K, n);
    return (SplusAJ (m[1], m[2]));
end;


SkewHermitianMatrixOctonionAlgebra := function (K, n)
    local m;
    m := MatrixOctonionAlgebra (K, n);
    return (SminusAJ (m[1], m[2]));
end;