/**********************************************************************/
/* Maxima file to check the derivation in "How to derive the field
   equations for three new kinds of gravitons"

   Version 1, Dec 3, 2024.
   Public domain, by David Parker, daveparker@pgu.org

   Command line used to run this file under Windows (the command lines
   for Macintosh and Linux should be similar):
     maxima -b gravitons_v1.mac

   File takes about 2.5 minutes to run on a laptop from 2020.

   Homepage for Maxima at:
     https://maxima.sourceforge.io/index.html
*/

display2d:false$

/**********************************************************************/
/* X (coordinate tensor), G (metric tensor), and inverse G (GInv). */

/* Maxima matrices are indexed starting from 1, but tensors are usually
   indexed starting from 0, so define matrix and tensor versions of
   X, G, and GInv.
*/

XMatrix: matrix( [ X0 ], [ X1 ], [ X2 ], [ X3 ] )$ /* Matrix. */
X(A) := XMatrix[A+1,1]$ /* Tensor. */

GMatrix: matrix( [ G00, G10, G20, G30 ] /* Matrix. */
               , [ G10, G11, G21, G31 ]
               , [ G20, G21, G22, G32 ]
               , [ G30, G31, G32, G33 ]
               )$
depends( [ G00, G10, G20, G30, G11, G21, G31, G22, G32, G33 ]
       , [ X0, X1, X2, X3 ]
       )$
G(A,B) := GMatrix[A+1,B+1]$ /* Tensor. */

GInvMatrix: matrix( [ GI00, GI10, GI20, GI30 ] /* Matrix. */
                  , [ GI10, GI11, GI21, GI31 ]
                  , [ GI20, GI21, GI22, GI32 ]
                  , [ GI30, GI31, GI32, GI33 ]
                  )$
depends( [ GI00, GI10, GI20, GI30, GI11, GI21, GI31, GI22, GI32, GI33 ]
       , [ X0, X1, X2, X3 ]
       )$
GInv(A,B) := GInvMatrix[A+1,B+1]$ /* Tensor. */

/**********************************************************************/
/* Equation (4), highest-order pure time (i.e. X0) derivatives. */

D: matrix( [ diff(G00,X0,1) ]
         , [ diff(G10,X0,1) ]
         , [ diff(G20,X0,1) ]
         , [ diff(G30,X0,1) ]
         , [ diff(G11,X0,2) ]
         , [ diff(G21,X0,2) ]
         , [ diff(G31,X0,2) ]
         , [ diff(G22,X0,2) ]
         , [ diff(G32,X0,2) ]
         , [ diff(G33,X0,2) ]
         )$

/**********************************************************************/
/* Equation (8): Christoffel symbol of the second kind: Gamma.

     ***** C
     *
     *
     *
     *
     * A B
*/

Gamma(A,B,C) := sum( GInv(C,D)*FirstKind(A,B,D), D, 0, 3 )$

/**********************************************************************/
/* Equation (10): Christoffel symbol of the first kind: [A B, C]. */

FirstKind(A,B,C) :=
  (1/2)*( diff(G(A,C),X(B)) + diff(G(B,C),X(A)) - diff(G(A,B),X(C)) )$

/**********************************************************************/
/* Equation (13): partial derivative of the inverse metric
   pd(g^{AB})/pd(x^C).
*/

PDGInvPDC( A, B, C ) :=
  - sum( sum( GInv(A,D)*GInv(B,E)*diff(G(D,E),X(C)), D, 0, 3 ), E, 0, 3 )$

/**********************************************************************/
/* Equation (7): Ricci tensor. */

RicciTensorElement(A,B) :=
    sum( diff(Gamma(A,B,C),X(C)), C, 0, 3 )
  - sum( diff(Gamma(A,C,C),X(B)), C, 0, 3 )
  + sum( sum( Gamma(A,B,D)*Gamma(C,D,C), C, 0, 3), D, 0, 3 )
  - sum( sum( Gamma(A,C,D)*Gamma(B,D,C), C, 0, 3), D, 0, 3 )$
RicciMatrix:
  genmatrix( lambda( [A,B], RicciTensorElement(A-1,B-1) ), 4, 4 )$
RicciTensor(A,B) := RicciMatrix[A+1,B+1]$

/* Replace derivatives of the inverse metric. */

for A: 0 thru 3 do
for B: 0 thru 3 do
for C: 0 thru 3 do
  RicciMatrix: subst( PDGInvPDC(A,B,C), diff(GInv(A,B),X(C)), RicciMatrix )$

/**********************************************************************/
/* Equation (6): Einstein tensor. */

EinsteinMatrix:
  genmatrix
  ( lambda
    ( [I,J]
    ,   RicciMatrix[I,J]
      - (1/2) * GMatrix[I,J]
      * sum( sum( GInvMatrix[K,L] * RicciMatrix[K,L], L, 1, 4), K, 1, 4)
    )
  , 4, 4
  )$
EinsteinTensor(A,B) := EinsteinMatrix[A+1,B+1]$

/**********************************************************************/
/* Equation (18): Q. */

Q( A, B, C, D ) :=
  factor
  ( sum
    ( sum
      ( GInv(E,F)
      * ( diff( FirstKind( A, B, F ), X(E) )
        - diff( FirstKind( A, E, F ), X(B) )
        )
      , F, D[1], D[2]
      )
    , E, C[1], C[2]
    )
  )$

/**********************************************************************/
/* Equation (19): P. */

P( A, B, C, D, E, F ) :=
  factor
  ( sum
    ( sum
      ( sum
        ( sum
          ( GInv(K,L)
          * GInv(M,N)
          * ( FirstKind( A, K, M )
            * FirstKind( B, L, N )
            - FirstKind( A, B, M )
            * FirstKind( K, L, N )
            )
          , N, F[1], F[2]
          )
        , M, E[1], E[2]
        )
      , L, D[1], D[2]
      )
    , K, C[1], C[2]
    )
  )$

/**********************************************************************/
/* Equation (20): check Ricci tensor equation using P and Q. */

for A: 0 thru 3 do
for B: 0 thru 3 do
  if( 0 # ratsimp( RicciTensor(A,B)
                 - ( Q(A,B,[0,3],[0,3]) + P(A,B,[0,3],[0,3],[0,3],[0,3]) )
                 )
    )
  then error( "ERROR: Ricci tensor equation (20) is incorrect." )$

"SUCCESS: Equation (20) for Ricci tensor checked."$

/**********************************************************************/
/* Equation (21): check Einstein tensor equation using P and Q. */

for A: 0 thru 3 do
for B: 0 thru 3 do
  if( 0 # ratsimp( EinsteinTensor(A,B)
                 - ( Q(A,B,[0,3],[0,3])
                   - (1/2)*G(A,B)
                     *sum( sum( GInv(C,D)*Q(C,D,[0,3],[0,3]), D,0,3), C,0,3)
                   + P(A,B,[0,3],[0,3],[0,3],[0,3])
                   - (1/2)*G(A,B)
                     *sum(sum( GInv(C,D)*P(C,D,[0,3],[0,3],[0,3],[0,3])
                             , D, 0, 3), C, 0, 3 )
                   )
                 )
    )
  then error( "ERROR: Einstein tensor equation (21) is incorrect." )$

"SUCCESS: Equation (21) for Einstein tensor checked."$

/**********************************************************************/
/* Theorem 1: check that HIP time (i.e. X0) derivatives are agree with
   equation (4).
*/

for I: 1 thru length(D) do
  if(   derivdegree(EinsteinMatrix,part(D[I,1],1),part(D[I,1],2))
      # part(D[I,1],3)
    )
  then error( "ERROR: HIP derivative", I, "incorrect" )$

"SUCCESS: Theorem 1 checked."$

/**********************************************************************/
/* Theorem 2: check linearity of Einstein matrix in HIP time
   derivatives.
*/

for I: 1 thru length(EinsteinMatrix) do
  for J: 1 thru I do
    block
    ( [Coeffs,Remainder,K]
    , Coeffs: zeromatrix(1,length(D))
    , for K: 1 thru length(D) do
        Coeffs[1,K]: coeff(EinsteinMatrix[I,J],D[K,1])
    , Remainder: EinsteinMatrix[I,J] - Coeffs . D
    , for K: 1 thru length(D) do
        if( not ( freeof(Coeffs,D[K,1]) and freeof(Remainder,D[K,1]) ) )
        then error( "ERROR: EinsteinMatrix[",I,",",J,"] not linear in D[",K,"]." )
    , if( 0 # ratsimp( EinsteinMatrix[I,J] - (Coeffs . D + Remainder) ) )
      then error( "ERROR: linearity test for Theorem 2 failed." )
    )$

"SUCCESS: Theorem 2 checked."$

/**********************************************************************/
/* Data and helper functions for checking the Appendices. */

/* Detect if a index is present in the P and Q arrays. */

IList: [0,I]$
JList: [0,J]$
KList: [0,K]$
LList: [0,L]$
MList: [0,M]$
NList: [0,N]$

/* Lexical sorting order for symbolic indices. */

Order[I]: 0$
Order[J]: 1$
Order[K]: 2$
Order[L]: 3$
Order[M]: 4$
Order[N]: 5$

/* Split ranges into 0 and [1,3]. */

Range: [ [0,0], [1,3] ]$

/* GInv(A,B), unordered indices. */
GIU(A,B) :=
  if( integerp(A) and integerp(B) ) then
    if (A<B) then GI(B,A) else GI(A,B)
  else
    if( ( A = 0 ) or ( ( B # 0 ) and ( Order[A] > Order[B] ) ) )
    then GI( B, A )
    else GI( A, B )$

/* First partial derivative of g_{ab}, unordered indices. */
PXU(A,B,C) :=
  if( integerp(A) and integerp(B) ) then
    if (A<B) then PX(B,A,C) else PX(A,B,C)
  else
    if( ( A = 0 ) or ( ( B # 0 ) and ( Order[A] > Order[B] ) ) )
    then PX( B, A, C )
    else PX( A, B, C )$

/* Second partial derivative of g_{ab}, unordered indices. */
PXXU(A,B,C,D) :=
  if( ( A = 0 ) or ( ( B # 0 ) and ( Order[A] > Order[B] ) ) )
  then
    if( ( D = 0 ) or ( ( C # 0 ) and ( Order[C] > Order[D] ) ) )
    then PXX( B, A, D, C )
    else PXX( B, A, C, D )
  else if( ( D = 0 ) or ( ( C # 0 ) and ( Order[C] > Order[D] ) ) )
    then PXX( A, B, D, C )
    else PXX( A, B, C, D )$

/* Christoffel symmbol of the 1st kind, unordered indices. */
C1U(A,B,C) := Half*(PXU(A,C,B)+PXU(B,C,A)-PXU(A,B,C))$

/* Partial derivative of C1U, unordered indices. */
PC1PXU(A,B,C,D) := Half*( PXXU(A,C,D,B) + PXXU(B,C,D,A) - PXXU(A,B,D,C))$

/**********************************************************************/
/* CheckExpanded -- in case the expressions in the Appendices for Q and P
   are not equal to the formulas for Q and P, expand the indices
   numerically to check for equality.  This is slower, so it is only
   tried if the symbolic check for equality fails.
*/

CheckExpanded( Expression, Formula, Prefix, ArgList, ArgNumList ) :=
  block
  ( [ExpandedExpression, ExpandedFormula, TempK, TempL, TempM, TempN]
  , ExpandedExpression:
      sum
      ( sum
        ( sum
          ( sum
            ( ev( Expression, K:TempK, L:TempL, M:TempM, N:TempN)
            , TempN, Range[ArgNumList[6]][1], Range[ArgNumList[6]][2]
            )
          , TempM, Range[ArgNumList[5]][1], Range[ArgNumList[5]][2]
          )
        , TempL, Range[ArgNumList[4]][1], Range[ArgNumList[4]][2]
        )
      , TempK, Range[ArgNumList[3]][1], Range[ArgNumList[3]][2]
      )
  , ExpandedFormula:
      sum
      ( sum
        ( sum
          ( sum
            ( ev( Formula, K:TempK, L:TempL, M:TempM, N:TempN)
            , TempN, Range[ArgNumList[6]][1], Range[ArgNumList[6]][2]
            )
          , TempM, Range[ArgNumList[5]][1], Range[ArgNumList[5]][2]
          )
        , TempL, Range[ArgNumList[4]][1], Range[ArgNumList[4]][2]
        )
      , TempK, Range[ArgNumList[3]][1], Range[ArgNumList[3]][2]
      )
  , ExpandedFormula: subst( [PX = PXU, GI=GIU], ExpandedFormula )
  , ExpandedFormula: ev( ExpandedFormula )
  , ExpandedFormula: ratsimp( ExpandedFormula )
  , return( ratsimp( ExpandedExpression - ExpandedFormula ) )
  )$

/**********************************************************************/
/* Check -- check the expressions in the Appendices for Q and P to see
   if they are equal to the formulas for Q and P.  They are first
   compared symbolically because this is relatively fast, but if that
   fails, call CheckExpanded to expand the indices numerically, which is
   slower.
*/

Check( Expression, Formula, Prefix, ArgList, ArgNumList ) :=
  if( 0 # ratsimp( Expression - Formula ) ) /* Try symbolic version first. */
  then
  ( /* If that didn't work, try expanding indices to integers. */
    if( 0 # CheckExpanded( Expression, Formula, Prefix, ArgList, ArgNumList ) ) then
      ( grind( "Expression: ", Expression )
      , grind( "Formula: ", Formula )
      , error( "ERROR in ", Prefix, ArgList )
      )
  )$

/*####################################################################*/
/* Appendix I. */

Q[0,0,0,0]: 0$
Q[0,0,0,L]: 0$
Q[0,0,K,0]: 0$
Q[0,0,K,L]: Half*GI(K,L)*(PXX(L,0,0,K)-PXX(0,0,K,L)-PXX(K,L,0,0)+PXX(K,0,0,L))$
Q[0,J,0,0]: 0$
Q[0,J,0,L]: Half*GI(L,0)*(PXX(J,L,0,0)-PXX(J,0,0,L)-PXX(L,0,0,J)+PXX(0,0,J,L))$
Q[0,J,K,0]: 0$
Q[0,J,K,L]: Half*GI(K,L)*(PXX(J,L,0,K)-PXX(J,0,K,L)-PXX(K,L,0,J)+PXX(K,0,J,L))$
Q[I,0,0,0]: 0$
Q[I,0,0,L]: 0$
Q[I,0,K,0]: Half*GI(K,0)*(PXX(0,0,I,K)-PXX(I,0,0,K)-PXX(K,0,0,I)+PXX(I,K,0,0))$
Q[I,0,K,L]: Half*GI(K,L)*(PXX(L,0,I,K)-PXX(I,0,K,L)-PXX(K,L,0,I)+PXX(I,K,0,L))$
Q[I,J,0,0]: Half*GI(0,0)*(PXX(J,0,0,I)-PXX(I,J,0,0)-PXX(0,0,I,J)+PXX(I,0,0,J))$
Q[I,J,0,L]: Half*GI(L,0)*(PXX(J,L,0,I)-PXX(I,J,0,L)-PXX(L,0,I,J)+PXX(I,0,J,L))$
Q[I,J,K,0]: Half*GI(K,0)*(PXX(J,0,I,K)-PXX(I,J,0,K)-PXX(K,0,I,J)+PXX(I,K,0,J))$
Q[I,J,K,L]: Half*GI(K,L)*(PXX(J,L,I,K)-PXX(I,J,K,L)-PXX(K,L,I,J)+PXX(I,K,J,L))$

QEquation(A,B,C,D) := GIU(C,D)*( PC1PXU(A,B,D,C) - PC1PXU(A,C,D,B) )$

for A: 1 thru 2 do
for B: 1 thru 2 do
for C: 1 thru 2 do
for D: 1 thru 2 do
  Check( Q[ IList[A], JList[B], KList[C], LList[D] ]
       , QEquation( IList[A], JList[B], KList[C], LList[D] )
       , 'Q
       , [IList[A],JList[B],KList[C],LList[D]]
       , [A,B,C,D,1,1]
       )$

"SUCCESS: Appendix I checked."$

/*####################################################################*/
/* Appendix II. */

P[0,0,0,0,0,0]: 0$
P[0,0,0,0,0,N]: 0$
P[0,0,0,0,M,0]: 0$
P[0,0,0,0,M,N]: 0$
P[0,0,0,L,0,0]: 0$
P[0,0,0,L,0,N]: 0$
P[0,0,0,L,M,0]: 0$
P[0,0,0,L,M,N]: 0$
P[0,0,K,0,0,0]: 0$
P[0,0,K,0,0,N]:Fourth*GI(K,0)*GI(N,0)
               *(PX(0,0,K)*(2*PX(N,0,0)-PX(0,0,N))
                -PX(0,0,0)*(PX(K,N,0)+PX(N,0,K)-PX(K,0,N))
                )$
P[0,0,K,0,M,0]:Fourth*GI(K,0)*GI(M,0)
               *((PX(M,0,K)+PX(K,M,0)-PX(K,0,M))*PX(0,0,0)
                -(2*PX(M,0,0)-PX(0,0,M))*PX(0,0,K)
                )$
P[0,0,K,0,M,N]:0$
P[0,0,K,L,0,0]:Fourth*GI(K,L)*GI(0,0)
               *(PX(0,0,K)*PX(0,0,L)
                -PX(0,0,0)*(PX(K,0,L)+PX(L,0,K)-PX(K,L,0))
                )$
P[0,0,K,L,0,N]:Fourth*GI(K,L)*GI(N,0)
               *(PX(0,0,K)*(PX(N,0,L)+PX(L,N,0)-PX(L,0,N))
                -PX(0,0,0)*(PX(K,N,L)+PX(L,N,K)-PX(K,L,N))
                )$
P[0,0,K,L,M,0]:Fourth*GI(K,L)*GI(M,0)
               *((PX(M,0,K)+PX(K,M,0)-PX(K,0,M))*PX(0,0,L)
                -(2*PX(M,0,0)-PX(0,0,M))*(PX(K,0,L)+PX(L,0,K)-PX(K,L,0))
                )$
P[0,0,K,L,M,N]:Fourth*GI(K,L)*GI(M,N)
               *((PX(M,0,K)+PX(K,M,0)-PX(K,0,M))*(PX(N,0,L)+PX(L,N,0)-PX(L,0,N))
                -(PX(K,N,L)+PX(L,N,K)-PX(K,L,N))*(2*PX(M,0,0)-PX(0,0,M))
                )$
P[0,J,0,0,0,0]:0$
P[0,J,0,0,0,N]:Fourth*GI(0,0)*GI(N,0)
               *(PX(0,0,0)*(PX(J,N,0)+PX(N,0,J)-PX(J,0,N))
                -PX(0,0,J)*(2*PX(N,0,0)-PX(0,0,N))
                )$
P[0,J,0,0,M,0]:Fourth*GI(0,0)*GI(M,0)
               *((2*PX(M,0,0)-PX(0,0,M))*PX(0,0,J)
                -(PX(M,0,J)+PX(J,M,0)-PX(J,0,M))*PX(0,0,0)
                )$
P[0,J,0,0,M,N]:0$
P[0,J,0,L,0,0]:Fourth*GI(L,0)*GI(0,0)
               *(PX(0,0,0)*(PX(J,0,L)+PX(L,0,J)-PX(J,L,0))
                -PX(0,0,J)*PX(0,0,L)
                )$
P[0,J,0,L,0,N]:Fourth*GI(L,0)*GI(N,0)
               *(PX(0,0,0)*(PX(J,N,L)+PX(L,N,J)-PX(J,L,N))
                -PX(0,0,J)*(PX(N,0,L)+PX(L,N,0)-PX(L,0,N))
                )$
P[0,J,0,L,M,0]:Fourth*GI(L,0)*GI(M,0)
               *((2*PX(M,0,0)-PX(0,0,M))*(PX(J,0,L)+PX(L,0,J)-PX(J,L,0))
                -(PX(M,0,J)+PX(J,M,0)-PX(J,0,M))*PX(0,0,L)
                )$
P[0,J,0,L,M,N]:Fourth*GI(L,0)*GI(M,N)
               *((2*PX(M,0,0)-PX(0,0,M))*(PX(J,N,L)+PX(L,N,J)-PX(J,L,N))
                -(PX(M,0,J)+PX(J,M,0)-PX(J,0,M))*(PX(N,0,L)+PX(L,N,0)-PX(L,0,N))
                )$
P[0,J,K,0,0,0]:0$
P[0,J,K,0,0,N]:Fourth*GI(K,0)*GI(N,0)
               *(PX(0,0,K)*(PX(J,N,0)+PX(N,0,J)-PX(J,0,N))
                -PX(0,0,J)*(PX(K,N,0)+PX(N,0,K)-PX(K,0,N))
                )$
P[0,J,K,0,M,0]:Fourth*GI(K,0)*GI(M,0)
               *((PX(M,0,K)+PX(K,M,0)-PX(K,0,M))*PX(0,0,J)
                -(PX(M,0,J)+PX(J,M,0)-PX(J,0,M))*PX(0,0,K)
                )$
P[0,J,K,0,M,N]:0$
P[0,J,K,L,0,0]:Fourth*GI(K,L)*GI(0,0)
               *(PX(0,0,K)*(PX(J,0,L)+PX(L,0,J)-PX(J,L,0))
                -PX(0,0,J)*(PX(K,0,L)+PX(L,0,K)-PX(K,L,0))
                )$
P[0,J,K,L,0,N]:Fourth*GI(K,L)*GI(N,0)
               *(PX(0,0,K)*(PX(J,N,L)+PX(L,N,J)-PX(J,L,N))
                -PX(0,0,J)*(PX(K,N,L)+PX(L,N,K)-PX(K,L,N))
                )$
P[0,J,K,L,M,0]:Fourth*GI(K,L)*GI(M,0)
               *((PX(M,0,K)+PX(K,M,0)-PX(K,0,M))*(PX(J,0,L)+PX(L,0,J)-PX(J,L,0))
                -(PX(M,0,J)+PX(J,M,0)-PX(J,0,M))*(PX(K,0,L)+PX(L,0,K)-PX(K,L,0))
                )$
P[0,J,K,L,M,N]:Fourth*GI(K,L)*GI(M,N)
               *((PX(M,0,K)+PX(K,M,0)-PX(K,0,M))*(PX(J,N,L)+PX(L,N,J)-PX(J,L,N))
                -(PX(M,0,J)+PX(J,M,0)-PX(J,0,M))*(PX(K,N,L)+PX(L,N,K)-PX(K,L,N))
                )$
P[I,0,0,0,0,0]:0$
P[I,0,0,0,0,N]:0$
P[I,0,0,0,M,0]:0$
P[I,0,0,0,M,N]:0$
P[I,0,0,L,0,0]:0$
P[I,0,0,L,0,N]:0$
P[I,0,0,L,M,0]:0$
P[I,0,0,L,M,N]:0$
P[I,0,K,0,0,0]:Fourth*GI(K,0)*GI(0,0)
               *((PX(I,0,K)+PX(K,0,I)-PX(I,K,0))*PX(0,0,0)
                -PX(0,0,I)*PX(0,0,K)
                )$
P[I,0,K,0,0,N]:Fourth*GI(K,0)*GI(N,0)
               *((PX(I,0,K)+PX(K,0,I)-PX(I,K,0))*(2*PX(N,0,0)-PX(0,0,N))
                -PX(0,0,I)*(PX(K,N,0)+PX(N,0,K)-PX(K,0,N))
                )$
P[I,0,K,0,M,0]:Fourth*GI(K,0)*GI(M,0)
               *((PX(I,M,K)+PX(K,M,I)-PX(I,K,M))*PX(0,0,0)
                -(PX(I,M,0)+PX(M,0,I)-PX(I,0,M))*PX(0,0,K)
                )$
P[I,0,K,0,M,N]:Fourth*GI(K,0)*GI(M,N)
               *((PX(I,M,K)+PX(K,M,I)-PX(I,K,M))*(2*PX(N,0,0)-PX(0,0,N))
                -(PX(I,M,0)+PX(M,0,I)-PX(I,0,M))*(PX(K,N,0)+PX(N,0,K)-PX(K,0,N))
                )$
P[I,0,K,L,0,0]:Fourth*GI(K,L)*GI(0,0)
               *((PX(I,0,K)+PX(K,0,I)-PX(I,K,0))*PX(0,0,L)
                -PX(0,0,I)*(PX(K,0,L)+PX(L,0,K)-PX(K,L,0))
                )$
P[I,0,K,L,0,N]:Fourth*GI(K,L)*GI(N,0)
               *((PX(I,0,K)+PX(K,0,I)-PX(I,K,0))*(PX(N,0,L)+PX(L,N,0)-PX(L,0,N))
                -PX(0,0,I)*(PX(K,N,L)+PX(L,N,K)-PX(K,L,N))
                )$
P[I,0,K,L,M,0]:Fourth*GI(K,L)*GI(M,0)
               *((PX(I,M,K)+PX(K,M,I)-PX(I,K,M))*PX(0,0,L)
                -(PX(I,M,0)+PX(M,0,I)-PX(I,0,M))*(PX(K,0,L)+PX(L,0,K)-PX(K,L,0))
                )$
P[I,0,K,L,M,N]:Fourth*GI(K,L)*GI(M,N)
               *((PX(I,M,K)+PX(K,M,I)-PX(I,K,M))*(PX(N,0,L)+PX(L,N,0)-PX(L,0,N))
                -(PX(I,M,0)+PX(M,0,I)-PX(I,0,M))*(PX(K,N,L)+PX(L,N,K)-PX(K,L,N))
                )$
P[I,J,0,0,0,0]:Fourth*GI(0,0)^2
               *(PX(0,0,I)*PX(0,0,J)
                -(PX(I,0,J)+PX(J,0,I)-PX(I,J,0))*PX(0,0,0)
                )$
P[I,J,0,0,0,N]:Fourth*GI(0,0)*GI(N,0)
               *(PX(0,0,I)*(PX(J,N,0)+PX(N,0,J)-PX(J,0,N))
                -(PX(I,0,J)+PX(J,0,I)-PX(I,J,0))*(2*PX(N,0,0)-PX(0,0,N))
                )$
P[I,J,0,0,M,0]:Fourth*GI(0,0)*GI(M,0)
               *((PX(I,M,0)+PX(M,0,I)-PX(I,0,M))*PX(0,0,J)
                -(PX(I,M,J)+PX(J,M,I)-PX(I,J,M))*PX(0,0,0)
                )$
P[I,J,0,0,M,N]:Fourth*GI(0,0)*GI(M,N)
               *((PX(I,M,0)+PX(M,0,I)-PX(I,0,M))*(PX(J,N,0)+PX(N,0,J)-PX(J,0,N))
                -(PX(I,M,J)+PX(J,M,I)-PX(I,J,M))*(2*PX(N,0,0)-PX(0,0,N))
                )$
P[I,J,0,L,0,0]:Fourth*GI(L,0)*GI(0,0)
               *(PX(0,0,I)*(PX(J,0,L)+PX(L,0,J)-PX(J,L,0))
                -(PX(I,0,J)+PX(J,0,I)-PX(I,J,0))*PX(0,0,L)
                )$
P[I,J,0,L,0,N]:Fourth*GI(L,0)*GI(N,0)
               *(PX(0,0,I)*(PX(J,N,L)+PX(L,N,J)-PX(J,L,N))
                -(PX(I,0,J)+PX(J,0,I)-PX(I,J,0))*(PX(N,0,L)+PX(L,N,0)-PX(L,0,N))
                )$
P[I,J,0,L,M,0]:Fourth*GI(L,0)*GI(M,0)
               *((PX(I,M,0)+PX(M,0,I)-PX(I,0,M))*(PX(J,0,L)+PX(L,0,J)-PX(J,L,0))
                -(PX(I,M,J)+PX(J,M,I)-PX(I,J,M))*PX(0,0,L)
                )$
P[I,J,0,L,M,N]:Fourth*GI(L,0)*GI(M,N)
               *((PX(I,M,0)+PX(M,0,I)-PX(I,0,M))*(PX(J,N,L)+PX(L,N,J)-PX(J,L,N))
                -(PX(I,M,J)+PX(J,M,I)-PX(I,J,M))*(PX(N,0,L)+PX(L,N,0)-PX(L,0,N))
                )$
P[I,J,K,0,0,0]:Fourth*GI(K,0)*GI(0,0)
               *((PX(I,0,K)+PX(K,0,I)-PX(I,K,0))*PX(0,0,J)
                -(PX(I,0,J)+PX(J,0,I)-PX(I,J,0))*PX(0,0,K)
                )$
P[I,J,K,0,0,N]:Fourth*GI(K,0)*GI(N,0)
               *((PX(I,0,K)+PX(K,0,I)-PX(I,K,0))*(PX(J,N,0)+PX(N,0,J)-PX(J,0,N))
                -(PX(I,0,J)+PX(J,0,I)-PX(I,J,0))*(PX(K,N,0)+PX(N,0,K)-PX(K,0,N))
                )$
P[I,J,K,0,M,0]:Fourth*GI(K,0)*GI(M,0)
               *((PX(I,M,K)+PX(K,M,I)-PX(I,K,M))*PX(0,0,J)
                -(PX(I,M,J)+PX(J,M,I)-PX(I,J,M))*PX(0,0,K)
                )$
P[I,J,K,0,M,N]:Fourth*GI(K,0)*GI(M,N)
               *((PX(I,M,K)+PX(K,M,I)-PX(I,K,M))*(PX(J,N,0)+PX(N,0,J)-PX(J,0,N))
                -(PX(I,M,J)+PX(J,M,I)-PX(I,J,M))*(PX(K,N,0)+PX(N,0,K)-PX(K,0,N))
                )$
P[I,J,K,L,0,0]:Fourth*GI(K,L)*GI(0,0)
               *((PX(I,0,K)+PX(K,0,I)-PX(I,K,0))*(PX(J,0,L)+PX(L,0,J)-PX(J,L,0))
                -(PX(I,0,J)+PX(J,0,I)-PX(I,J,0))*(PX(K,0,L)+PX(L,0,K)-PX(K,L,0))
                )$
P[I,J,K,L,0,N]:Fourth*GI(K,L)*GI(N,0)
               *((PX(I,0,K)+PX(K,0,I)-PX(I,K,0))*(PX(J,N,L)+PX(L,N,J)-PX(J,L,N))
                -(PX(I,0,J)+PX(J,0,I)-PX(I,J,0))*(PX(K,N,L)+PX(L,N,K)-PX(K,L,N))
                )$
P[I,J,K,L,M,0]:Fourth*GI(K,L)*GI(M,0)
               *((PX(I,M,K)+PX(K,M,I)-PX(I,K,M))*(PX(J,0,L)+PX(L,0,J)-PX(J,L,0))
                -(PX(I,M,J)+PX(J,M,I)-PX(I,J,M))*(PX(K,0,L)+PX(L,0,K)-PX(K,L,0))
                )$
P[I,J,K,L,M,N]:Fourth*GI(K,L)*GI(M,N)
               *((PX(I,M,K)+PX(K,M,I)-PX(I,K,M))*(PX(J,N,L)+PX(L,N,J)-PX(J,L,N))
                -(PX(I,M,J)+PX(J,M,I)-PX(I,J,M))*(PX(K,N,L)+PX(L,N,K)-PX(K,L,N))
                )$

PEquation( A, B, C, D, E, F ) :=
  GIU(C,D)*GIU(E,F)*( C1U(A,C,E)*C1U(B,D,F) - C1U(A,B,E)*C1U(C,D,F) )$

for A: 1 thru 2 do
for B: 1 thru 2 do
for C: 1 thru 2 do
for D: 1 thru 2 do
for E: 1 thru 2 do
for F: 1 thru 2 do
  block
  ( [Temp]
  , Temp: subst( Half^2, Fourth, P[ IList[A], JList[B], KList[C], LList[D], MList[E], NList[F] ] )
  , Check( Temp
         , PEquation( IList[A], JList[B], KList[C], LList[D], MList[E], NList[F] )
         , 'P
         , [IList[A],JList[B],KList[C],LList[D],MList[E],NList[F]]
         , [A,B,C,D,E,F]
         )
  )$

"SUCCESS: Appendix II checked."$

/*####################################################################*/

"SUCCESS: All checks succeeded."$