/* Tensor package. This code intends to simplify tensorial expressions.
 */

/* functions internal to tensors */
RuleBase("Delta",{ind1,ind2});
RuleBase("TList",{head,tail});
RuleBase("TSum",{indices,body});
RuleBase("TD",{ind});
RuleBase("X",{ind});

/* And the simplificaiton rules for X, addition, subtraction
   and multiplication */
10 # (TD(_i)X(_j)) <-- Delta(i,j);
10 # (TD(_i) ( (_f) + (_g) ) ) <-- (TD(i)f) + (TD(i)g);
10 # (TD(_i) ( (_f) - (_g) ) ) <-- (TD(i)f) - (TD(i)g);
10 # (TD(_i) (      - (_g) ) ) <--          -  TD(i)g;
10 # (TD(_i) ( (_f) * (_g) ) ) <-- (TD(i)f)*g + f*(TD(i)g);
10 # (TD(_i) ( (_f) ^ (n_IsPositiveInteger) ) ) <-- n*(TD(i)f)*f^(n-1);
10 # (TD(_i)Delta(_j,_k)) <-- 0;
10 # (TD(_i)f_IsNumber) <-- 0;

/* The only TSum summation simplification: summing over no indices
   means no summation. */
10 # (TSum({})(_body)) <-- body;

/* Explicit summation when Ndim is defined. This summation will
   be invoked when using TExplicitSum. */
20 # (TSum(_indices)(_body))_(IsInteger(Ndim)) <--
    LocalSymbols(index,i,sum)
    [
      Local(index,i,sum);
      index:=indices[1];
      sum:=0;
      MacroLocal(index);
      For(i:=1,i<=Ndim,i++)
      [
        MacroSet(index,i);
        sum:=sum+Eval(TSum(Tail(indices))body);
      ];
      sum;
    ];

/* TExplicitSum sets the dimension of the space under consideration,
   so summation can proceed */
(TExplicitSum(Ndim_IsInteger)(_body)) <-- Eval(body);

/* Move the delta factors to the front, so they can be simplified
   away. It uses ApplyDelta to move a factor either to the front
   or to the back of the list. Input is a list of factors, as
   returned by Flatten(expressions,"*")
   */
MoveDeltas(_list) <--
[
  Local(result,i,nr);
  result:={};
  nr:=Length(list);
  For(i:=1,i<=nr,i++)
  [
    ApplyDelta(result,list[i]);
  ];
  result;
];


10  # ApplyDelta(_result,Delta(_i,_j)) <--
    DestructiveInsert(result,1,Delta(i,j));
20 # ApplyDelta(_result,(_x) ^ (n_IsInteger))_(n>0) <--
    [
      Local(i);
      For(i:=1,i<=n,i++)
      [
        ApplyDelta(result,x);
      ];
    ];
100 # ApplyDelta(_result,_term) <--
    DestructiveAppend(result,term);
    

/* TSimplify : expand brackets, and send the expression of addition
   of terms to TSimplifyAux */
TSimplify(TSum(_indices)(_f)) <--
[
  TSimplifyAux(TSum(indices)ExpandBrackets(f));
];


/* TSimplifyAux : simplify each term independently */
10 # TSimplifyAux(TSum(_indices)((_f) + (_g))) <--
     TSimplifyAux(TSum(FlatCopy(indices))(f)) +
     TSimplifyAux(TSum(FlatCopy(indices))(g));
10 # TSimplifyAux(TSum(_indices)((_f) - (_g))) <--
     TSimplifyAux(TSum(FlatCopy(indices))(f)) -
     TSimplifyAux(TSum(FlatCopy(indices))(g));
10 # TSimplifyAux(TSum(_indices)(   - (_g))) <--
                                    -  TSimplifyAux(TSum(indices)(g));

40 # TSimplifyAux(TSum(_indices)_body) <--
[
  Local(flat);

  /* Convert expressions of the form (a*b*c) to {a,b,c} */
  flat:=Flatten(body,"*");

  /* Move the deltas to the front. */
  flat:=MoveDeltas(flat);

  /* Simplify the deltas away (removing the required indices) */
  flat:=TSumRest(flat);

  /* Determine if there are indices the summand still depends on */
  Local(varlist,independ,nrdims);
  varlist:=VarList(flat);
  independ:=Intersection(indices,varlist);
  nrdims:=Length(indices)-Length(independ);

  /* Return result, still summing over the indices not removed by deltas */
  Ndim^nrdims*TSum(independ)flat;
];


/* Terminating condition for the tensorial simplification */

10 # TSumSimplify(TList(Delta(_ind,_ind),_list))_Contains(indices,ind) <--  

[
  /* Remove the index from the list of indices to sum over, since
     it is now implicitly summed over by simplifying the delta */
  DestructiveDelete(indices,Find(indices,ind));

/* Return result simplified for this delta */
  Ndim*TSumRest(list);
];

11 # TSumSimplify(TList(Delta(_ind1,_ind2),_list))_
     Contains(indices,ind2) <--
[
  /* Remove the index from the list of indices to sum over, since
     it is now implicitly summed over by simplifying the delta */
  DestructiveDelete(indices,Find(indices,ind2));

  /* Return result simplified for this delta */
  TSumRest( Subst(ind2,ind1)list );
];
11 # TSumSimplify(TList(Delta(_ind1,_ind2),_list))_
     Contains(indices,ind1) <--
[
  /* Remove the index from the list of indices to sum over, since
     it is now implicitly summed over by simplifying the delta */
  DestructiveDelete(indices,Find(indices,ind1));

  /* Return result simplified for this delta */
  TSumRest( Subst(ind1,ind2)list );
];



1010 # TSumSimplify(TList(_term,_list)) <--
[
  term*TSumRest(list);
];


10 # TSumRest({}) <-- 1;
20 # TSumRest(_list) <--
[
   TSumSimplify(TList(Head(list),Tail(list)));
];



UnFence("TSumSimplify",1); 
UnFence("TSumRest",1); 
UnFence("TSum",2);