homkil.mac

/* number of columns of a matrix */
numcol (m) :=
    if length (m) = 0 then 0
        else length (m [1]);

/* flatten matrix into a column vector */
flatmatc (m) :=
    block (
        [a],
        a : matrix (),
        for i : 1 thru numcol (m) do
            a : addrow (a, col (m, i)),
        a
    );

/* flatten list of matrices into a long-row matrix */
flatmatlsr (l) :=
    block (
        [a],
        a : matrix (),
        for i : 1 thru length (l) do
            a : addcol (a, flatmatc (l [i])),
        a
    );

/* flatten list of matrices into a long-column matrix */
flatmatlsc (l) :=
    block (
        [a],
        a : matrix (),
        for i : 1 thru length (l) do
            a : addrow (a, l [i]),
        a
    );

/* flatten matrix of matrices into a long-column matrix */
flatblockmatc (m) :=
    block (
        [a, ri],
        a : matrix (),
        for i : 1 thru length (m) do (
            ri : m [i],
            for j : 1 thru length (ri) do
                a : addrow (a, ri [j])
        ),
        a
    );

/* multiply a matrix by a list of matrices */
matmulls (m, l) :=
    block (
        [a],
        a : [],
        for i :1 thru length (l) do
            a : endcons (m . l [i], a),
        a
    );

/* multiply 2 block matrices */
blockmatmul (m1, m2) :=
    map (
        lambda ([r1],
            apply ("+",
                map (
                    lambda ([e1, r2],
                        map (lambda ([e2], e1 . e2), r2)
                    ),
                    r1, args (m2)
                )
            )
        ),
        m1
    );

/* linear combination given by a column vector */
collincomb (c, b) :=
    block (
        [a],
        a : 0,
        for i : 1 thru length (b) do
            a : a + c [i] [1] * b [i],
        a
    );

/* least-square solution of m . x = y */
lssolve (m) :=
    block (
        [mct],
        mct : conjugate (transpose (m)),
        ((mct . m)^^-1) . mct
    );

/* adjoint representation of a matrix Lie algebra */
adjrep2 (l1, l2) :=
    block (
        [n1, n2, l2ss, a, aa, x, y, c],
        n1 : length (l1),
        n2 : length (l2),
        l2ss : lssolve (flatmatlsr (l2)),
        a : [],
        for i : 1 thru n1 do (
            x : l1 [i],
            aa : matrix (),
            for j : 1 thru n2 do (
                y : l2 [j],
                c : x . y - y . x,
                aa : addcol (aa, l2ss . flatmatc (c))
            ),
            a : endcons (aa, a)
        ),
        a
    );

adjrep (l) :=
    adjrep2 (l, l);

/* "curvature" of a Lie algebra map, i.e. matrices
   r (x) . r (y) - r (y) . r (x) - r (ad (x) (y)) */
lacurv (r, ad) :=
    block (
        [n, a, aa, rx, ry, adx, c1, c2],
        n : length (r),
        a : matrix (),
        for i : 1 thru n do (
            rx : r [i],
            adx : ad [i],
            aa : [],
            for j : 1 thru n do (
                ry : r [j],
                c1 : rx . ry - ry . rx,
                c2 : collincomb (col (adx, j), r),
                aa : endcons (c1 - c2, aa)
            ),
            a : addrow (a, aa)
        ),
        a
    );

/* check Lie algebra structure */
checkla (l, ad) :=
    block (
        [c],
        c : lacurv (l, ad),
        if zeromatrixp (c) then return (true),
        print (box ("Not a Lie algebra, obstruction:")),
        print (c),
        false
    );

/* Killing form of a Lie algebra */
lakil (l) :=
    apply ('matrix,
        map (
            lambda ([a],
                map (
                    lambda ([b],
                        mattrace (a . b)
                    ),
                    l
                )
            ),
            l
        )
    );

/* standard 2nd order Casimir element */
casimir2 (r, ad) :=
    block (
        [k],
        k : lakil (ad),
        first (first (
            blockmatmul (
                matrix (r),
                (k^^-1) . covect (r)
            )
        ))
    );

/* direct sum of 2 Lie algebras */
lasum (l1, l2) :=
    block (
        [n1, n2, z1, z2],
        n1 : length (first (l1)),
        n2 : length (first (l2)),
        z1 : zeromatrix (n1, n1),
        z2 : zeromatrix (n2, n2),
        append (
            map (lambda ([m1], diag ([m1, z2])), l1),
            map (lambda ([m2], diag ([z1, m2])), l2)
        )
    );

/* split list to fist n elements and the rest */
nsplitls (l, n) :=
    block (
        [m],
        m : length (l),
        if m <= n
            then [l, []]
            else [firstn (l, n), lastn (l, m - n)]
    );

/* list of n zero matrices n x n */
zeromatls (n) :=
    makelist (zeromatrix (n, n), n);

/* elementary basis vector, or zero if k > n */
evect (k, n) :=
    if k <= n
        then ematrix (n, 1, 1, k, 1)
        else zeromatrix (n, 1);

/* affine LA matrix from evect and a matrix */
affmat (k, m) :=
    block (
        [nr, nc],
        nr : length (m),
        nc : numcol (m),
        addrow (
            zeromatrix (1, 1 + nc),
            addcol (evect (k, nr), m)
        )
    );

/* list of affine LA matrices */
affmatls (l) :=
    block (
        [a],
        a : [],
        for i : 1 thru length (l) do
            a : endcons (affmat (i, l [i]), a),
        a
    );

/* decompose matrix to 2 blocks horizontally */
blockmat2h (m, nc) :=
    block (
        [a, b],
        a : matrix (),
        b : matrix (),
        for i : 1 thru numcol (m) do
            if i <= nc
                then a : addcol (a, col (m, i))
                else b : addcol (b, col (m, i)),
        matrix ([a, b])
    );

/* decompose matrix to 2 blocks vertically */
blockmat2v (nr, m) :=
    block (
        [a, b],
        a : matrix (),
        b : matrix (),
        for i : 1 thru length (m) do
            if i <= nr
                then a : addrow (a, row (m, i))
                else b : addrow (b, row (m, i)),
        matrix ([a], [b])
    );

/* decompose matrix to four blocks */
blockmat4 (nr, m, nc) :=
    block (
        [mv],
        mv : blockmat2v (nr, m),
        addrow (
            blockmat2h (mv [1] [1], nc),
            blockmat2h (mv [2] [1], nc)
        )
    );

/* decompose list of matrices to four blocks */
blockmat4ls (nr, l, nc) :=
    block (
        [a11, a12, a21, a22, b],
        a11 : [],
        a12 : [],
        a21 : [],
        a22 : [],
        for i : 1 thru length (l) do (
            b : blockmat4 (nr, l [i], nc),
            a11 : endcons (b [1] [1], a11),
            a12 : endcons (b [1] [2], a12),
            a21 : endcons (b [2] [1], a21),
            a22 : endcons (b [2] [2], a22)
        ),
        matrix ([a11, a12], [a21, a22])
    );

/* decompose matrix of matrices to four blocks */
blockmat4mat (nr, m, nc) :=
    block (
        [a11, a12, a21, a22, bl],
        a11 : matrix (),
        a12 : matrix (),
        a21 : matrix (),
        a22 : matrix (),
        for i : 1 thru length (m) do (
            bl : blockmat4ls (nr, m [i], nc),
            a11 : addrow (a11, bl [1] [1]),
            a12 : addrow (a12, bl [1] [2]),
            a21 : addrow (a21, bl [2] [1]),
            a22 : addrow (a22, bl [2] [2])
        ),
        matrix ([a11, a12], [a21, a22])
    );

/* Levi-Civita connection from torsion */
lcconn1 (t, g, i, j) :=
    block (
        [gi, tj],
        gi : g [i],
        tj : apply (addcol, t [j]),
        transpose (gi . tj . (g^^-1))
    );

lcconn (t, g) :=
    block (
        [n, a, aa, c, tij, gtij, gtji],
        n : length (t),
        a : [],
        for i : 1 thru n do (
            aa : matrix (),
            for j : 1 thru n do (
                tij : t [i] [j],
                gtij : lcconn1 (t, g, i, j),
                gtji : lcconn1 (t, g, j, i),
                c : (gtij + gtji - tij) / 2,
                aa : addcol (aa, c)
            ),
            a : endcons (aa, a)
        ),
        a
    );

/* check affine structure (= reductive pair) */
checkaffcurv (blcurv) :=
    block (
        [c0],
        c0 : matrix (
            [zerofor (blcurv [1] [1]),
                blcurv [1] [2]],
            [blcurv [2] [1],
                blcurv [2] [2]]
        ),
        if zeromatrixp (c0) then return (true),
        print (box ("Not an affine structure, obstruction:")),
        print (c0),
        false
    );

/* apply a representation on a metric (bilinear form) */
repmet (r, g) :=
    block (
        [n, a, c, gc],
        n : length (g),
        a : [],
        for i : 1 thru length (r) do (
            c : blockmat4 (n, r [i], n) [1] [1],
            gc : transpose (c) . g + g . c,
            a : endcons (gc, a)
        ),
        a    
    );

/* check invariant metric (bilinear form) */
checkinvmet (r, g) :=
    block (
        [rg],
        rg : repmet (r, g),
        if zeromatrixp (rg) then return (true),
        print (box ("Not an invariant metric, obstruction:")),
        print (rg),
        false
    );

/* check Levi-Civita connection (= torsion free) */
checklctors (tors) :=
    block (
        if zeromatrixp (tors) then return (true),
        print (box ("Not a Levi-Civita connection, torsion:")),
        print (tors),
        false
    );

/* prepare Riemannian homogeneous space from [l, g];
   l = basis of a LA split to V + h
   g = metric on V

   returns [g, ad, lc, curv];
   ad = adjoint representation of l
   lc = Levi-Civita connection
   curv = curvature of lc */
riemhom (l, g) :=
    block (
        [n, ad, adspl, hc, aff, blcurv, acurv, tors, lc, curv],
        n : length (g),
        ad : adjrep (l),
        if not checkla (l, ad) then return (false),
        adspl : nsplitls (ad, n),
        hc : blockmat4ls (n, adspl [2], n) [1] [1],
        if not checkinvmet (hc, g) then return (false),
        aff : affmatls (append (zeromatls (n), hc)),
        blcurv : blockmat4 (n, lacurv (aff, ad), n),
        if not checkaffcurv (blcurv) then return (false),
        acurv : blockmat4mat (1, blcurv [1] [1], 1),
        tors : acurv [2] [1],
        lc : append (lcconn (tors, g), hc),
        aff : affmatls (lc),
        blcurv : blockmat4 (n, lacurv (aff, ad), n),
        acurv : blockmat4mat (1, blcurv [1] [1], 1),
        tors : acurv [2] [1],
        if not checklctors (tors) then return (false),
        curv : acurv [2] [2],
        [g, ad, lc, curv]
    );

/* Ricci tensor from curvature */
ricci (curv) :=
    block (
        [n],
        n : length (curv),
        apply ('matrix,
            makelist (
                makelist (
                    sum (curv[k,i][k,j], k, 1, n),
                    j, n
                ),
                i, n
            )
        )
    );

/* solve overdetermined system by integrability conditions */
lasolve (r, ad) :=
    block (
        [m, m1, o, l, s, ls, rr],
        print ("Computing LA curvature"),
        m : flatblockmatc (lacurv (r, ad)),
        print ("Elimination pass", 0),
        m : gcdelim (m),
        o : 1,
        do (
            print ("Multiplication pass", o),
            m1 : addrow (m, flatmatlsc (matmulls (m, r))),
            print ("Elimination pass", o),
            m1 : gcdelim (m1),
            if length (m) = length (m1) then return ('done),
            m : m1,
            o : o + 1
        ),
        print ("Solving"),
        /* empty matrix does not know number of columns */
        if emptyp (m) then
            [s, l] : [ident (length (first (r))), 1]
        else
            [s, l] : gcdsolve (m),
        print ("Restricting representation"),
        if numcol (s) > 0 then (
            ls : lssolve (s),
            rr : map (lambda ([e], ratsimp (ls . e . s)), r)
        ) else rr : map (lambda ([e], 0), r),
        [o, l, s, rr]
    );

/* dual, i.e. negative transpose of a matrix */
dualrep (m) :=
    - transpose (m);

/* dual LA representation */
dualrepls (l) :=
    map (dualrep, l);

/* dual of matrix of matrices */
dualrepmat (m) :=
    matrixmap (dualrep, m);

/* add first index to exterior multi-index */
addextidx (i1, l) :=
    block (
        [a],
        a : [i1],
        for i : 1 thru length (l) do
            a : endcons (l [i] + i1, a),
        a
    );

/* remove first index from exterior multi-index */
remextidx (l) :=
    block (
        [i1, a],
        i1 : l [1],
        a : [],
        for i : 2 thru length (l) do
            a : endcons (l [i] - i1, a),
        a
    );

/* replace index in exterior multi-index */
repextidx (l, r, i) :=
    block (
        [l1, p, i1, j, j1, s],
        l1 : copylist (l),
        i1 : l1 [r],
        s : 1,
        if i1 > i then (
            j : r,
            do (
                j1 : j,
                if j1 = 1 then return ('done),
                j : j - 1,
                i1 : l1 [j],
                if i1 < i then return ('done),
                if i1 = i then (
                    s : 0,
                    return ('done)
                ),
                l1 [j1] : i1,
                s : -s
            ),
            l1 [j1] : i
        ) else if i1 < i then (
            p : length (l),
            j : r,
            do (
                j1 : j,
                if j1 = p then return ('done),
                j : j + 1,
                i1 : l1 [j],
                if i1 > i then return ('done),
                if i1 = i then (
                    s : 0,
                    return ('done)
                ),
                l1 [j1] : i1,
                s : -s
            ),
            l1 [j1] : i
        ),
        [l1, s]
    );

/* merge index into multi-index by value */
extpidx (l, i) :=
    block (
        [p, j, j1, i1, a, s],
        p : length (l),
        s : 1,
        j : 0,
        do (
            j1 : j,
            if j1 = p then return ('done),
            j : j + 1,
            i1 : l [j],
            if i1 > i then return ('done),
            if i1 = i then (
                s : 0,
                return ('done)
            ),
            s : -s
        ),
        a : append (
            firstn (l, j1),
            [i],
            lastn (l, p - j1)
        ),
        [a, s]
    );

/* find and remove index from multi-index by value */
intpidx (l, i) :=
    block (
        [p, a, s],
        p : length (l),
        a : [],
        s : 0,
        for j : 1 thru p do (
            if l [j] = i then (
                a : append (
                    firstn (l, j - 1),
                    lastn (l, p - j)
                ),
                s : (-1)^ (j - 1)
            )
        ),
        [a, s]
    );

/* split flat index to exterior multi-index */
splitextidx (k, n, p) :=
    block (
        [m, l, k1],
        if p = 0 then return ([]),
        m : binomial (n, p) - k,
        l : p,
        do (
            k1 : binomial (l, p) - m,
            if k1 > 0 then return ('done),
            l : l + 1
        ),
        addextidx (n - l + 1, splitextidx (k1, l - 1, p - 1))
    );

/* flatten exterior multi-index */
flatextidx (l, n) :=
    block (
        [p, n1, k1],
        p : length (l),
        if p = 0 then return (1),
        n1 : n - l [1],
        k1 : binomial (n, p) - binomial (n1 + 1, p),
        k1 + flatextidx (remextidx (l), n1)
    );

/* exterior power of a matrix (in LA sense) */
extrep (m, p) :=
    block (
        [n, np, a, mi, mi1, s, k1, v, x],
        n : length (m),
        np : binomial (n, p),
        a : zeromatrix (np, np),
        for k : 1 thru np do (
            mi : splitextidx (k, n, p),
            for r : 1 thru p do (
                v : col (m, mi [r]),
                for l : 1 thru n do (
                    x : v [l] [1],
                    mi1 : repextidx (mi, r, l),
                    s : mi1 [2],
                    if s # 0 then (
                        k1 : flatextidx (mi1 [1], n),
                        a [k1] [k] : a [k1] [k] + s * x
                    )
                )
            )
        ),
        a
    );

/* exterior power of a list of matrices */
extrepls (l, p) :=
    block (
        [a, b],
        a : [],
        for i : 1 thru length (l) do (
            b : extrep (l [i], p),
            a : endcons (b, a)
        ),
        a
    );

/* exterior power of a matrix of matrices */
extrepmat (m, p) :=
    block (
        [a, b],
        a : matrix (),
        for i : 1 thru length (m) do (
            b : extrepls (m [i], p),
            a : addrow (a, b)
        ),
        a
    );

/* matrix of exterior product */
extp (v, p) :=
    block (
        [n, p1, np, np1, a, mi, mi1, s, k1, x],
        n : length (v),
        p1 : p + 1,
        np : binomial (n, p),
        np1 : binomial (n, p1),
        a : zeromatrix (np1, np),
        for k : 1 thru np do (
            mi : splitextidx (k, n, p),
            for l : 1 thru n do (
                x : v [l] [1],
                mi1 : extpidx (mi, l),
                s : mi1 [2],
                if s # 0 then (
                    k1 : flatextidx (mi1 [1], n),
                    a [k1] [k] : a [k1] [k] + s * x
                )
            )
        ),
        a
    );

/* matrix of interior product */
intp (v, p) :=
    block (
        [n, p1, np, np1, a, mi, mi1, s, k1, x],
        n : length (v),
        p1 : p - 1,
        np : binomial (n, p),
        np1 : binomial (n, p1),
        a : zeromatrix (np1, np),
        for k : 1 thru np do (
            mi : splitextidx (k, n, p),
            for l : 1 thru n do (
                x : v [l] [1],
                mi1 : intpidx (mi, l),
                s : mi1 [2],
                if s # 0 then (
                    k1 : flatextidx (mi1 [1], n),
                    a [k1] [k] : a [k1] [k] + s * x
                )
            )
        ),
        a
    );

/* partially skew-symmetrized action of curvature */
curvext1 (curvps, v, p, ids) :=
    block (
        [n, vc],
        n : length (curvps),
        vc : transpose (v) . curvps,
        apply ("+", makelist (
            block (
                [ei],
                ei : evect (i, n),
                tenprod (extp (ei, p), ids) . (vc . ei)
            ),
            i, n
        ))
    );

/* totally skew-symmetrized action of curvature */
curvext2 (curvps, p, ids) :=
    block (
        [n, p1],
        n : length (curvps),
        p1 : p + 1,
        apply ("+", makelist (
            block (
                [ei, eic],
                ei : evect (i, n),
                eic : transpose (ei) . curvps,
                tenprod (extp (ei, p1), ids)
                . apply ("+", makelist (
                    block (
                        [ej],
                        ej : evect (j, n),
                        tenprod (extp (ej, p), ids) . (eic . ej)
                    ),
                    j, n
                ))
            ),
            i, n
        ))
    );

/* representation on Killing forms */
kfrep (conn, curv, p) :=
    block (
        [m, n, p1, ids,
        connd, connp, connp1, curvd, curvp,
        a, b, b11, b12, b21, b22, ei],
        print ("Preparing Killing representation"),
        m : length (conn),
        n : length (curv),
        p1 : p + 1,
        ids : ident (1),
        connd : dualrepls (conn),
        connp : extrepls (connd, p),
        connp1 : extrepls (connd, p1),
        curvd : dualrepmat (curv),
        curvp : extrepmat (curvd, p),
        a : [],
        for i : 1 thru m do (
            ei : evect (i, n),
            b22 : connp [i],
            b11 : connp1 [i],
            b21 : - intp (ei, p1),
            b12 : - 1 / p * curvext1 (curvp, ei, p, ids),
            b : addrow (
                addcol (b11, b12),
                addcol (b21, b22)
            ),
            a : endcons (b, a)
        ),
        a
    );

/* solve Killing forms by integrability conditions */
kfsolve (rh, p) :=
    block (
        [g, ad, conn, curv],
        [g, ad, conn, curv] : firstn (rh, 4),
        lasolve (kfrep (conn, curv, p), ad)
    );

/* check complex structure */
checkcplx (J) :=
    block (
        [IJ],
        IJ : ident (length (J)) + J . J,
        if zeromatrixp (IJ) then return (true),
        print (box ("Not a complex structure, obstruction:")),
        print (IJ),
        false
    );

/* contact distribution and its complement */
contactsplit (J) :=
    block (
        [J2, IJ2, d, z],
        J2 : J . J,
        IJ2 : ident (length (J)) + J2,
        d : transpose (gcdelim (transpose (-J2))),
        z : transpose (gcdelim (transpose (IJ2))),
        [IJ2, d, z]
    );

/* check contact structure */
checkcontact (J, IJ2, d, z) :=
    block (
        [n, dn, zn, Jdz],
        [n, dn, zn] : [length (J), numcol (d), numcol (z)],
        if zn # 1 or zn + dn # n then (
            print (box ("Wrong dimensions of contact structure:")),
            print (n, dn, zn),
            return (false)
        ),
        Jdz : addcol (IJ2 . d, J . z),
        if zeromatrixp (Jdz) then return (true),
        print (box ("Not a contact structure, obstruction:")),
        print (Jdz),
        false
    );

/* complementary maximal isotropic subspaces */
isotrop (J, d) :=
    block (
        [i, dp, dn],
        i : %i * ident (length (J)),
        dp : transpose (gcdelim (transpose ((J + i) . d / 2))),
        dn : transpose (gcdelim (transpose ((J - i) . d / 2))),
        [dp, dn]
    );

/* check Hermitian or contact metric */
checkhermit (g, J) :=
    block (
        [gJ],
        gJ : transpose (J) . g + g . J,
        if zeromatrixp (gJ) then return (true),
        print (box ("Not a Hermitian or contact metric, obstruction:")),
        print (gJ),
        false
    );

/* count one bits in integer */
bit_countones (i, m) :=
    block (
        [c],
        c : 0,
        for j : 0 while j < m do
            if bit_onep (i, j) then c : c + 1,
        c
    );

/* Clifford algebra in standard basis */
stdcliffpos (i, m, bm) :=
    block (
        [bi],
        i : i - 1,
        bi : bit_lsh (1, i),
        genmatrix (
            lambda ([j, k],
                k : k - 1,
                if bit_xor (j - 1, k, bi) = 0
                        and bit_and (k, bi) = 0 then
                    - (-1) ^ bit_countones (k, i)
                else 0
            ),
            bm, bm
        )
    );

stdcliffneg (i, m, bm) :=
    block (
        [bi],
        i : i - 1,
        bi : bit_lsh (1, i),
        genmatrix (
            lambda ([j, k],
                k : k - 1,
                if bit_xor (j - 1, k, bi) = 0
                        and bit_and (k, bi) # 0 then
                    2 * (-1) ^ bit_countones (k, i)
                else 0
            ),
            bm, bm
        )
    );

stdcliffzero (m, bm) :=
    genmatrix (
        lambda ([j, k],
            if j = k then
                %i * (-1) ^ bit_countones (k - 1, m)
            else 0
        ),
        bm, bm
    );

stdcliffeven (n) :=
    block (
        [m, bm],
        m : n / 2,
        bm : bit_lsh (1, m),
        genmatrix (
            lambda ([j, k],
                if k <= m then stdcliffpos (k, m, bm)
                else stdcliffneg (k - m, m, bm)
            ),
            1, n
        )
    );

stdcliffodd (n) :=
    block (
        [m2, m],
        m2 : n - 1,
        m : m2 / 2,
        bm : bit_lsh (1, m),
        genmatrix (
            lambda ([j, k],
                if k <= m then stdcliffpos (k, m, bm)
                elseif k <= m2 then stdcliffneg (k - m, m, bm)
                else stdcliffzero (m, bm)
            ),
            1, n
        )
    );

/* Clifford algebra for general metric */
cliffeven (n, g, J) :=
    block (
        [dp, dn, dntg, cp, cn],
        if not checkhermit (g, J) then return (false),
        if not checkcplx (J) then return (false),
        [dp, dn] : isotrop (J, ident (n)),
        dntg : transpose (dn) . g,
        cp : ((dntg . dp) ^^ -1) . dntg,
        cn : transpose (dp) . g,
        stdcliffeven (n) . addrow (cp, cn)
    );

cliffodd (n, g, J) :=
    block (
        [IJ2, d, z, ztg, dp, dn, dntg, cp, cn, c0],
        if not checkhermit (g, J) then return (false),
        [IJ2, d, z] : contactsplit (J),
        if not checkcontact (J, IJ2, d, z) then return (false),
        [dp, dn] : isotrop (J, d),
        dntg : transpose (dn) . g,
        cp : ((dntg . dp) ^^ -1) . dntg,
        cn : transpose (dp) . g,
        ztg : transpose (z) . g,
        c0 : ztg / sqrt (ztg . z),
        stdcliffodd (n) . addrow (cp, cn, c0)
    );

cliff (g, J) :=
    block (
        [n],
        n : length (g),
        if evenp (n) then cliffeven (n, g, J)
        else cliffodd (n, g, J)
    );

/* prepare spin homogeneous space from [rh, J];
   r = [g, ad, conn, curv] = Riemannian homogeneous space
   J = complex or contact structure for constructing spinors

   returns [g, ad, lc, curv, c];
   c = Clifford algebra */
spinhom (rh, J) :=
    append (rh, [cliff (first (rh), J)]);

/* spin representation of a matrix (in LA sense) */
spinrep (m, g, c) :=
    blockmatmul (c, - m . (g ^^ -1) . transpose (c) / 4) [1, 1];

spinrepls (l, g, c) :=
    block (
        [gc],
        gc : - (g ^^ -1) . transpose (c) / 4,
        map (lambda ([m], blockmatmul (c, m . gc) [1, 1]), l)
    );

spinrepmat (b, g, c) :=
    block (
        [gc],
        gc : - (g ^^ -1) . transpose (c) / 4,
        matrixmap (lambda ([m], blockmatmul (c, m . gc) [1, 1]), b)
    );

/* representation on Killing spinors */
ksrep (conn, g, c, a) :=
    block (
        [m, n, conns],
        print ("Preparing Killing representation"),
        m : length (conn),
        n : length (curv),
        conns : spinrepls (conn, g, c),
        makelist (
            conns [i]
            - if i <= n then a * c [1, i] else 0,
            i, m
        )
    );

/* solve Killing spinors by integrability conditions */
kssolve (sh, a) :=
    block (
        [g, ad, conn, curv, c],
        [g, ad, conn, curv, c] : firstn (sh, 5),
        lasolve (ksrep (conn, g, c, a), ad)
    );

/* tensor product of two matrices */
tenprod (m1, m2) :=
    apply (matrix,
        apply (append,
            map (
                lambda ([r1],
                    map (
                        lambda ([r2],
                            apply (append,
                                map (
                                    lambda ([e1],
                                        map (
                                            lambda ([e2],
                                                e1 * e2
                                            ),
                                            r2
                                        )
                                    ),
                                    r1
                                )
                            )
                        ),
                        args (m2)
                    )
                ),
                args (m1)
            )
        )
    );

/* tensor product representation (in LA sense) */
tenprodrep (m1, m2) :=
    tenprod (m1, ident (length (m2)))
    + tenprod (ident (length (m1)), m2);

tenprodrepls (l1, l2) :=
    map (tenprodrep, l1, l2);

tenprodrepmat (b1, b2) :=
    map (tenprodrepls, b1, b2);

/* skew-symmetrized square of Killing homomorphism */
kilh2ext (A) :=
    block (
        [A2],
        A2 : blockmatmul (transpose (A), A),
        A2 - transpose (A2)
    );

/* shift Kiling homomorphism up to order p + 1 */
kilhup1 (A, v, p, ids) :=
    block (
        [n, p1, Av, iv],
        n : numcol (A),
        p1 : p + 1,
        Av : A . v,
        iv : tenprod (intp (v, p1), ids),
        1 / p * apply ("+", makelist (
            block (
                [ei],
                ei : evect (i, n),
                tenprod (extp (ei, p), ids) . (
                    Av . tenprod (intp (ei, p1), ids)
                    - (A . ei) . iv
                )
            ),
            i, n
        ))
    );

kilhup1ext (A, p, ids) :=
    block (
        [n, p1],
        n : numcol (A),
        p1 : p + 1,
        apply ("+", makelist (
            block (
                [ei],
                ei : evect (i, n),
                tenprod (extp (ei, p1), ids)
                . kilhup1 (A, ei, p, ids)
            ),
            i, n
        ))
    );

kilhup (A, v, p, ids) :=
    block (
        [n, p1, p2, Av],
        n : numcol (A),
        p1 : p + 1,
        p2 : p + 2,
        kilhup1 (A, v, p, ids)
        - (
            if p2 > n then 0
            else 1 / 2 / p1 * (
                tenprod (intp (v, p2), ids)
                . kilhup1ext (A, p, ids)
            )
        )
    );

/* representation on Killing spinor-valued forms */
ksfrep (conn, curv, g, c, p, A) :=
    block (
        [m, n, p1, p2, ids,
        connd, connp, connp1, conns, connps, connp1s,
        curvd, curvp, curvs, curvps, khcurvps, khce2,
        a, b, b11, b12, b21, b22, ei],
        print ("Preparing Killing representation"),
        m : length (conn),
        n : length (curv),
        p1 : p + 1,
        p2 : p + 2,
        ids : ident (length (c [1, 1])),
        connd : dualrepls (conn),
        connp : extrepls (connd, p),
        connp1 : extrepls (connd, p1),
        conns : spinrepls (conn, g, c),
        connps : tenprodrepls (connp, conns),
        connp1s : tenprodrepls (connp1, conns),
        curvd : dualrepmat (curv),
        curvp : extrepmat (curvd, p),
        curvs : spinrepmat (curv, g, c),
        curvps : tenprodrepmat (curvp, curvs),
        khcurvps : curvps + kilh2ext (A),
        khce2 : 1 / 2 / p1 * curvext2 (khcurvps, p, ids),
        a : [],
        for i : 1 thru m do (
            ei : evect (i, n),
            b22 : connps [i] - A . ei,
            b11 : connp1s [i] - kilhup (A, ei, p, ids),
            b21 : - tenprod (intp (ei, p1), ids),
            b12 : - 1 / p * (
                curvext1 (khcurvps, ei, p, ids)
                - (
                    if p2 > n then 0
                    else tenprod (intp (ei, p2), ids) . khce2
                )
            ),
            b : addrow (
                addcol (b11, b12),
                addcol (b21, b22)
            ),
            a : endcons (b, a)
        ),
        a
    );

/* solve Killing spinor-valued forms by integrability conditions */
ksfsolve (sh, p, A) :=
    block (
        [g, ad, conn, curv, c],
        [g, ad, conn, curv, c] : firstn (sh, 5),
        lasolve (ksfrep (conn, curv, g, c, p, A), ad)
    );

/* standard Killing homomorphism */
stdkilh (c, p) :=
    block (
        [n, idp],
        n : numcol (c),
        idp : ident (binomial (n, p)),
        matrix (map (
            lambda ([e], tenprod (idp, e)),
            first (c)
        ))
    );

/* solve standard Killing spinor-valued forms */
stdksfsolve (sh, p, a) :=
    block (
        [g, ad, conn, curv, c, A],
        [g, ad, conn, curv, c] : firstn (sh, 5),
        A : a * stdkilh (c, p),
        ksfsolve (sh, p, A)
    );

/* raise degree by Clifford multiplication */
clraise (c, p) :=
    block (
        [n],
        n : numcol (c),
        apply ("+", makelist (
            block (
                [ei],
                ei : evect (i, n),
                tenprod (extp (ei, p), c [1, i])
            ),
            i, n
        ))
    );

/* lower degree by Clifford multiplication */
cllower (c, g, p) :=
    block (
        [cg, n],
        cg : c . (g^^-1),
        n : numcol (cg),
        apply ("+", makelist (
            block (
                [ei],
                ei : evect (i, n),
                tenprod (intp (ei, p), c [1, i])
            ),
            i, n
        ))
    );

/* elementary skew symmetric matrix */
esomat (n, x, i, j) :=
    ematrix (n, n, x, i, j) - ematrix (n, n, x, j, i);

/* orthogonal LA */
sola (n) :=
    apply (append,
        makelist (
            makelist (
                esomat (n, 1, i, j),
                i, j + 1, n
            ),
            j, n
        )
    );

/* standard sphere */
sphere (n) :=
    block (
        [l, g, J],
        l : sola (n + 1),
        g : ident (n),
        J : stdJ (n),
        ratsimp (spinhom (riemhom (l, g), J))
    );

/* euclidean LA */
eucla (n) :=
    block (
        [n1],
        n1 : n + 1,
        apply (append,
            cons (
                makelist (
                    ematrix (n1, n1, 1, i, 1),
                    i, 2, n1
                ),
                makelist (
                    makelist (
                        esomat (n1, 1, i, j),
                        i, j + 1, n1
                    ),
                    j, 2, n1
                )
            )
        )
    );

/* euclidean space */
euclid (n) :=
    block (
        [l, g, J],
        l : eucla (n),
        g : ident (n),
        J : stdJ (n),
        ratsimp (spinhom (riemhom (l, g), J))
    );

/* elementary symmetric matrix */
esymat (n, x, i, j) :=
    ematrix (n, n, x, i, j) + ematrix (n, n, x, j, i);

/* Lorentzian LA */
lorla (n) :=
    apply (append,
        cons (
            makelist (
                esymat (n, 1, i, 1),
                i, 2, n
            ),
            makelist (
                makelist (
                    esomat (n, 1, i, j),
                    i, j + 1, n
                ),
                j, 2, n
            )
        )
    );

/* hyperbolic space */
hyperb (n) :=
    block (
        [l, g, J],
        l : lorla (n + 1),
        g : ident (n),
        J : stdJ (n),
        ratsimp (spinhom (riemhom (l, g), J))
    );

/* elementary unitary matrix */
eumat (n, x, i, j) :=
    block (
        [em],
        em : ematrix (n, n, x, i, j),
        if i = j then em : em / 2,
        em - conjugate (transpose (em))
    );

/* unitary LA */
ula (n) :=
    block (
        [a],
        a : [],
        for i : 1 thru n do (
            a : endcons (eumat (n, %i, i, i), a),
            for j : i + 1 thru n do (
                a : endcons (eumat (n, 1, j, i), a),
                a : endcons (eumat (n, %i, j, i), a)
            )
        ),
        a
    );

/* special unitary LA */
sula (n) :=
    block (
        [a],
        a : [],
        for i : 1 thru n - 1 do (
            a : endcons (
                eumat (n, %i, i, i) - eumat (n, %i, n, n),
                a
            )
        ),
        for i : 1 thru n do (
            for j : i + 1 thru n do (
                a : endcons (eumat (n, 1, j, i), a),
                a : endcons (eumat (n, %i, j, i), a)
            )
        ),
        a
    );

/* symplectic LA */
spla1 () :=
    [
        matrix ([-%i, 0], [0, %i]),
        matrix ([0, -1], [1, 0]),
        matrix ([0, %i], [%i, 0])
    ];

spla2 () :=
    block (
        [sp1],
        sp1 : spla1 (),
        append (
            [
                kronecker_product (matrix ([0, -1], [1, 0]), ident (2))
            ],
            map (
                lambda ([i%],
                    kronecker_product (matrix ([0, 1], [1, 0]), i%)
                ),
                sp1
            ),
            map (
                lambda ([i%],
                    kronecker_product (matrix ([0, 0], [0, 1]), i%)
                ),
                sp1
            ),
            map (
                lambda ([i%],
                    kronecker_product (matrix ([1, 0], [0, 0]), i%)
                ),
                sp1
            )
        )
    );

/* standard complex or contact structure */
stdJ2 : matrix ([0, 1], [-1, 0]);

stdJ (n) :=
    block (
        [m, J],
        m : floor (n / 2),
        J : tenprod (stdJ2, ident (m)),
        if oddp (n) then J : diag ([0, J]),
        J
    );

/* Berger sphere */
bergs (m, s) :=
    block (
        [l, n, g, J],
        l : ula (m),
        n : 2 * m - 1,
        g : ident (n) + ematrix (n, n, s - 1, 1, 1),
        J : stdJ (n),
        ratsimp (spinhom (riemhom (l, g), J))
    );

/* quaternionic sphere */
quats2 (r, s) :=
    block (
        [sp2, l, g, J],
        sp2 : spla2 (),
        l : append (
            [sp2[5], sp2[6], sp2[1], sp2[2], sp2[7], sp2[3], sp2[4]],
            rest (sp2, 7)
        ),
        g : diag ([r, r, s, s, r, s, s]),
        J : stdJ (7),
        ratsimp (spinhom (riemhom (l, g), J))
    );

/* Aloff-Wallach space W^{1,1} */
w11 (r, s) :=
    block (
        [u3, l, g, J],
        u3 : ula (3),
        l : [
            u3[1] - u3[6],
            u3[2], u3[4], u3[5],
            -u3[3], -u3[7], -u3[8],
            -u3[1] - u3[6] + 2 * u3[9]
        ],
        g : diag ([r, r, s, s, r, s, s]),
        J : stdJ (7),
        ratsimp (spinhom (riemhom (l, g), J))
    );