The 4-dimensional traces in FORM are rather sophisticated. The manual 
explains what FORM does exactly but sometimes people do not appreciate what 
the consequences are. Let us first have a look at something seemingly 
unrelated. If one is in 4 dimensions a totally antisymmetric object with 5 
indices must be zero. Hence e_(m1,m2,m3,m4,m5) = 0. And therefore also 
e_(m1,m2,m3,m4,m5)*e_(n1,n2,n3,n4,n5) must be zero (the m and n are 
4-dimensional indices). This means that if we contract those two 
Levi-Civita tensors, we get 120 terms, each the product of 5 Kronecker 
delta's and together those terms add up to zero. Similarly when the m and n 
are 4-vectors we get 120 terms, each with 5 dotproducts and the sum is 
zero. That means that there are nontrivial relations between those 
apparently random objects.

The basis relatie is the Schouten identiteit:
  e_(m1,m2,m3,m4)*d_(m5,m6) =
		+e_(m5,m2,m3,m4)*d_(m1,m6)
		+e_(m1,m5,m3,m4)*d_(m2,m6)
		+e_(m1,m2,m5,m4)*d_(m3,m6)
		+e_(m1,m2,m3,m5)*d_(m4,m6)
This identity is however amazingly difficult to apply in practise, because 
it is usually not clear what to take for the part outside the Levi-Civita 
tensor. This holds in particular when substitutions like p1.p1 = m1^2 have 
been made (or even p1.p1 = 0).

Let us look now at a simple example: We take the trace of 6 gamma matrices 
together with a gamma5. Cut out the next program and run it.
%-----------------------------g6.frm--------------------------------------
I	m1,m2,m3,m4,m5,m6,m7,m8,m9,m10;
I	n1,n2,n3,n4,n5,n6;
CF	g5,d;
.global
*
*	The expression F is the expression that will give the FORM answer.
*	The expression G is what most other programs produce when this trace
*	is taken.
*
L	F = g_(1,5_,m1,m2,m3,m4,m5,m6);
L	G = g5(m1,m2,m3,m4,m5,m6);
id	g5(?a) = distrib_(-1,2,d,g5,?a);
id	d(m1?,m2?) = d_(m1,m2);
id	g5(m1?,m2?,m3?,m4?) = g_(1,5_,m1,m2,m3,m4);
Trace4,1;
Print +f +s;
.sort
*
*	Here we substitute values for the indices. Effectively this gives all
*	combinations although it may not look like it.
*	We see that the expressions are eventually always the same.
*
Multiply replace_(m1,1,m2,2,m3,3,m4,4);
Print +f +s;
.sort
Drop;
#do m5 = 1,4
#do m6 = 1,4
L	F`m5'`m6' = (F - G)*replace_(m5,`m5',m6,`m6');
#enddo
#enddo
Print +f +s;
.store
*
*	Restart with F and G
*
L	F = g_(1,5_,m1,m2,m3,m4,m5,m6);
L	G = g5(m1,m2,m3,m4,m5,m6);
id	g5(?a) = distrib_(-1,2,d,g5,?a);
id	d(m1?,m2?) = d_(m1,m2);
id	g5(m1?,m2?,m3?,m4?) = g_(1,5_,m1,m2,m3,m4);
Trace4,1;
*
*	Now we apply 6 times the Schouten identity to force the expressions into
*	the same format. Eventually we get the same formula.
*
#do i = 1,6
id	e_(n1?,n2?,n3?,n4?)*d_(m`i',n5?) =
		+e_(m`i',n2,n3,n4)*d_(n1,n5)
		+e_(n1,m`i',n3,n4)*d_(n2,n5)
		+e_(n1,n2,m`i',n4)*d_(n3,n5)
		+e_(n1,n2,n3,m`i')*d_(n4,n5);
#enddo
Print +f +s;
.store
*
*	A demonstration how this can run out of hand when there are more gamma's
*	First with 8
*
L	F = g_(1,5_,m1,m2,m3,m4,m5,m6,m7,m8);
L	G = g5(m1,m2,m3,m4,m5,m6,m7,m8);
repeat id g5(m1?,m2?,m3?,m4?,?a) = distrib_(-1,2,d,g5,m1,m2,m3,m4,?a);
id	d(m1?,m2?) = d_(m1,m2);
id	g5(m1?,m2?,m3?,m4?) = g_(1,5_,m1,m2,m3,m4);
Trace4,1;
Print +f +s;
.store
*
*	A demonstration how this can run out of hand when there are more gamma's
*	Now with 10.
*
L	F = g_(1,5_,m1,m2,m3,m4,m5,m6,m7,m8,m9,m10);
L	G = g5(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10);
repeat id g5(m1?,m2?,m3?,m4?,?a) = distrib_(-1,2,d,g5,m1,m2,m3,m4,?a);
id	d(m1?,m2?) = d_(m1,m2);
id	g5(m1?,m2?,m3?,m4?) = g_(1,5_,m1,m2,m3,m4);
Trace4,1;
.end
%-----------------------------g6.frm--------------------------------------
Now let us look at a realistic example. We look at the reaction
  e-e+ -> gamma* -> tau-tau+ -> nutau d ubar nubartau mu+ numu
This is a 2->6 reaction but there is only one Feynman diagram.
We program it twice, once with the contract option in the trace (default) 
and once without:
%-----------------------------tau.frm-------------------------------------
#-
S	me,mt,mnt,md,mu,mnm,mm;
I	m1,m2,m3,m4,m5,m6;
V	pe,pp,pa,pb,pnt,pntb,pd,pub,pnm,pm;
V	p1,p2,p3,p4;

Local F1 = (g_(1,pp)-me)*g_(1,m1)*(g_(1,pe)+me)*g_(1,m2)
	*(g_(2,pnt)+mnt)*g_(2,m3,7_)*(g_(2,pa)+mt)*g_(2,m1)
	*(g_(2,pb)+mt)*g_(2,m5,7_)*(g_(2,pntb)+mnt)
	*g_(2,m6,7_)*(g_(2,pb)+mt)*g_(2,m2)*(g_(2,pa)+mt)*g_(2,m4,7_)
		*(g_(3,pd)+md)*g_(3,m3,7_)*(g_(3,pub)+mu)*g_(3,m4,7_)
		*(g_(4,pnm)+mnm)*g_(4,m5,7_)*(g_(4,pm)+mm)*g_(4,m6,7_)
		;
Local F2 = (g_(1,pp)-me)*g_(1,m1)*(g_(1,pe)+me)*g_(1,m2)
	*(g_(2,pnt)+mnt)*g_(2,m3,7_)*(g_(2,pa)+mt)*g_(2,m1)
	*(g_(2,pb)+mt)*g_(2,m5,7_)*(g_(2,pntb)+mnt)
	*g_(2,m6,7_)*(g_(2,pb)+mt)*g_(2,m2)*(g_(2,pa)+mt)*g_(2,m4,7_)
		*(g_(3,pd)+md)*g_(3,m3,7_)*(g_(3,pub)+mu)*g_(3,m4,7_)
		*(g_(4,pnm)+mnm)*g_(4,m5,7_)*(g_(4,pm)+mm)*g_(4,m6,7_)
		;
*Local G = (g_(1,pp)-me)*g_(1,m1)*(g_(1,pe)+me)*g_(1,m2)
*	*(g_(2,pnt)+mnt)*g_(2,m3)*(gi_(2)-g_(2,5_))*(g_(2,pa)+mt)*g_(2,m1)
*	*(g_(2,pb)+mt)*g_(2,m5)*(gi_(2)-g_(2,5_))*(g_(2,pntb)+mnt)
*	*g_(2,m6)*(gi_(2)-g_(2,5_))*(g_(2,pb)+mt)*g_(2,m2)
*	*(g_(2,pa)+mt)*g_(2,m4)*(gi_(2)-g_(2,5_))
*		*(g_(3,pd)+md)*g_(3,m3)*(gi_(3)-g_(3,5_))
*		*(g_(3,pub)+mu)*g_(3,m4)*(gi_(3)-g_(3,5_))
*		*(g_(4,pnm)+mnm)*g_(4,m5)*(gi_(4)-g_(4,5_))
*		*(g_(4,pm)+mm)*g_(4,m6)*(gi_(4)-g_(4,5_))
*		;
if ( expression(F1,G) );
	Trace4,4;
	Trace4,3;
	Trace4,2;
	Trace4,1;
else;
	Trace4,nocontract,4;
	Trace4,nocontract,3;
	Trace4,nocontract,2;
	Trace4,nocontract,1;
	Contract,0;
endif;
Print +f +s;
.sort
*Drop G;
id	p1?.p2? = p1(0)*p2(0)-p1(1)*p2(1)-p1(2)*p2(2)-p1(3)*p2(3);
id	e_(p1?,p2?,p3?,p4?) = -e_(p1,p2,p3,p4)*e_(0,1,2,3);
Contract;
Print +f +s;
.sort
Drop;
L	FF = F1-F2;
Print +f +s;
.end
%-----------------------------tau.frm-------------------------------------
When you run this program, you will see that the first version (F1) is done 
very fast and gives only 27 terms. The second version does not use contract 
and produces a few thousand terms, many with a Levi-Civita tensor. Because 
we work here in a theory without CP violation, the answer should be time 
reversal invariant and hence the Levi-Civita tensor terms should add up to 
zero. To prove this with the Schouten identity is rather hard. But also the 
terms without Levi-Civita terms are different. The second expression has 
very many of them. Here the first remark about the contraction of 
Levi-Civita tensors with 5 indices applies. So how can we see that the two 
expressions are actually identical in 4 dimensions? This is done by NOT 
making any substitutions of dotproducts yet but writing everything in 
vector components. We obtain then two equally long expressions and when we 
subtract them the answer is zero. Hence these expressions must be identical 
numerically. Had we however first made substitutions on them things might 
become very difficult.

If you want to have some fun, uncomment the expression G in the program 
tau.frm and run the trace again. Now g7_ -> 1-g5_ and there are 256 times 
as many terms. The answer is the same, but there is quite some work to be 
done. A program like REDUCE has to take its traces in this way. 
Additionally, last time I used REDUCE, it did not know the Chisholm 
identities and after a few hours we had to kill the program (we had not put 
the neutrino and quark masses in it to make it simpler while the above 
program has them and hence does more work).

The problem with other programs is that their trace algorithms are 
different from those in FORM. This means that in 4 dimensions they will 
almost always give a seemingly different result. But that does not mean 
that either one is wrong. It is just that usually FORM is more efficient 
and that it is difficult to prove that the results are identical.

For n-dimensional traces the story is different. There we cannot apply 
identities as the Chisholm identity and gamma5 does not exist technically. 
In REDUCE it will still take the trace with the algorithm that we used for 
expression G in the first program. In FORM the philosophy is that you have 
to provide your own gamma5 algorithm, which may be the one we used for the 
G expression or maybe something else. For example, in the calculation of 
three loop anomalous dimensions and coefficient functions we eliminated the 
gamma5 with the relation:
    g5_(1) = e_(m1,m2,m3,m4)*g_(1,m1,m2,m3,m4)/24
after which we applied a finite renormalization. The Levi-Civita tensor 
will combine for instance with another Levi-Civita tensor from a projector 
and the two are contracted with their indices being in n dimensions. Taking 
the algorithm for expression G in the first program would not work well 
here. Hence, be on your guard with gamma5 in n dimensions. Do not use built 
in algorithms of any program without checking and knowing what you are 
doing.