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

   Public domain, by David Parker, daveparker@pgu.org

   Version 1, Dec 3, 2024.
   Version 2, Dec 23, 2024: updated some of the comments.
   Version 3, Jan 10, 2025: updated some of the comments, changed order
                            of the derivatives in D.

   Command line used to run this file under Windows (the command lines
   for Macintosh and Linux should be similar):
     maxima -b gravitons_v3.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. */

/**********************************************************************/
/* \label{eq:HighestOrderDerivatives}, 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(G22,X0,2) ]
         , [ diff(G33,X0,2) ]
         , [ diff(G21,X0,2) ]
         , [ diff(G31,X0,2) ]
         , [ diff(G32,X0,2) ]
         )$

/**********************************************************************/
/* label{eq:ChristoffelSecondKind}: 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 )$

/**********************************************************************/
/* \label{eq:InverseMetricDerivatives}: 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 )$

/**********************************************************************/
/* \label{eq:ChristoffelFirstKind}: 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)) )$

/**********************************************************************/
/* \label{eq:RicciTensorChristoffelSecondKind}: 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 )$

/**********************************************************************/
/* \label{eq:EinsteinTensorREliminated}: 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]$

/**********************************************************************/
/* \label{eq:QDefinition}: 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]
    )
  )$

/**********************************************************************/
/* \label{eq:PDefinition}: 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]
    )
  )$

/**********************************************************************/
/* \label{eq:ExplicitRicciTensor}: 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: eq:ExplicitRicciTensor is incorrect." )$

"SUCCESS: eq:ExplicitRicciTensor checked."$

/**********************************************************************/
/* \label{eq:ExplicitEinsteinTensor}: 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: eq:ExplicitEinsteinTensor is incorrect." )$

"SUCCESS: eq:ExplicitEinsteinTensor checked."$

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

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."$