https://github.com/vermaseren/form/pull/524#issuecomment-2124965359

test.frm

#define SIZE "5"
*#include- include/color.h
#define RANK "6"
*
*	The file with all the declarations for running color traces.
*	The variable RANK tells how many Adjoint invariants of how many
*	non-adjoint invariants can occur simultaneously. For 14 vertices
*	this number is maximally 4 (4 loops in a non-adjoint representation)
*
*	Note: One needs version 3 of FORM to run this.
*	Version 2 (or 2.3) makes no chance. The program would need very
*	extensive reworking and it would become far less efficient.
*	All internal variables start with the string cOl to avoid potential 
*	conflicts with namings in the calling program.
*
*	The external variables are:
*	Tensor T;
*	Tensor f(antisymmetric);
*	Symbol NA;
*	Symbol NR;
*	Symbol I2R;
*	Symbol cA;
*	Symbol cR;
*	Symbol [cR-cA/2];
*	All generic d contractions like d33, d763 etc. which are CFunctions.
*
*	File by J.Vermaseren, 24-may-1997
*
*--#[ Declarations :
*
Symbol NA,NR,I2R;
Dimension NA;
AutoDeclare Index cOli,cOlj,cOlk,cOln;
AutoDeclare Symbol cOlI;
AutoDeclare Vector cOlp,cOlq;
AutoDeclare Symbol cOlx,cOly,cOlc;
AutoDeclare Tensor cOld;
AutoDeclare Tensor cOldr(symmetric),cOlda(symmetric);
Vector cOlpR1,...,cOlpR`RANK',cOlpA1,...,cOlpA`RANK';
Tensor cOldR1,...,cOldR`RANK',cOldA1,...,cOldA`RANK';
Tensor,<cOldr1(symmetric)>,...,<cOldr`RANK'(symmetric)>
		,<cOlda1(symmetric)>,...,<cOlda`RANK'(symmetric)>;
Symbol cOln, cA,cR,[cR-cA/2];
Tensor cOlfp,cOlff(cyclic),T,cOlTr(cyclic),cOlTt(cyclic),cOlf3;
NTensor cOlTN,cOlTM;
CFunction cOlTT,cOlnum,cOldddff,cOldff554,cOld,cOldd,cOlddd,cOlorig,cOlacc;
CFunction cOlE,cOlE0,cOlEa,cOlEb,cOlEc,cOldexp;
CFunction d33,d44,d55,d433,d66,d633,d543,d444,d3333
	,d77,d743,d653,d644,d554,d5333,d4433a,d4433b,d4433c
	,d88,d853,d844,d763,d754,d7333,d664,d655,d6433a,d6433b,d6433c
	,d5533a,d5533b,d5533c,d5533d,d5443a,d5443b,d5443c,d5443d
	,d4444a,d4444b,d43333a,d43333b;
Tensor f(antisymmetric),cOlf1(antisymmetric);
Set cOlpAs:cOlpA1,...,cOlpA`RANK';
Set cOlpRs:cOlpR1,...,cOlpR`RANK';
Set cOlpAR:cOlpA1,...,cOlpA`RANK',cOlpR1,...,cOlpR`RANK';
Set cOldar:cOlda1,...,cOlda`RANK',cOldr1,...,cOldr`RANK';
Set cOldas:cOlda1,...,cOlda`RANK';
Set cOliii:cOli1,...,cOli30;
Set cOljjj:cOlj1,...,cOlj30;
ProperCount ON;
*
*--#] Declarations : 
*--#[ SORT :
*
*	Next Procedure handles the .sort. This way one can put the program
*	in a checking mode in which there is printing at all .sort instructions
*	by just adding a single statement.
*
#procedure SORT(text)
*Print +f +s;
.sort:`text';
#endprocedure
*
*--#] SORT : 
*--#[ adjoint :
*
#procedure adjoint
*
*	Procedure to deal with gluonic loops (adjoint representation)
*	In this case there are special shortcuts for odd loops.
*	Also even loops have special savings.
*	Use the declarations of the file cfactor.h
*	Do not use indices cOlk,cOlk1,cOlk2,cOlk3,cOlk4,cOlk5 when calling this routine!
*
*	Routine by J.Vermaseren 24-may-1997
*
#define adj "0"
repeat;
	if ( count(cOlff,1) == 0 );
		ReplaceLoop,f,a=3,l=all,outfun=cOlff;
		id	cOlff(cOli1?,cOli2?) = -cA*d_(cOli1,cOli2);
		id	cOlff(cOli1?,cOli2?,cOli3?) = cA/2*f(cOli1,cOli2,cOli3);
	endif;
endrepeat;
if ( count(cOlff,1) ) redefine adj "1";
renumber;
#call SORT(adjoint-1)
#if `adj' > 0
#do iStageB = 1,1
	id,once,cOlff(cOli1?,cOli2?,cOli3?,cOli4?,cOli5?,cOli6?,cOli7?) = 1/2*(
		+cOlff(cOlk1,cOli6,cOli5,cOli4,cOli3,cOli2)*f(cOli1,cOli7,cOlk1)
		+cOlff(cOlk1,cOli5,cOli4,cOli3,cOli2,cOli7)*f(cOli1,cOli6,cOlk1)
		+cOlff(cOli1,cOli5,cOli4,cOli3,cOli2,cOlk1)*f(cOli7,cOli6,cOlk1)
		+cOlff(cOli1,cOlk1,cOli3,cOli2,cOli6,cOli7)*f(cOli5,cOli4,cOlk1)
		+cOlff(cOli1,cOli4,cOlk1,cOli2,cOli6,cOli7)*f(cOli5,cOli3,cOlk1)
		+cOlff(cOli1,cOli4,cOli3,cOlk1,cOli6,cOli7)*f(cOli5,cOli2,cOlk1)
		+cOlff(cOli1,cOlk1,cOli2,cOli5,cOli6,cOli7)*f(cOli4,cOli3,cOlk1)
		+cOlff(cOli1,cOli3,cOlk1,cOli5,cOli6,cOli7)*f(cOli4,cOli2,cOlk1)
		+cOlff(cOli1,cOlk1,cOli4,cOli5,cOli6,cOli7)*f(cOli3,cOli2,cOlk1)
	);
#call SORT(adjoint-2-a)
	id,once,cOlff(cOli1?,cOli2?,cOli3?,cOli4?,cOli5?,cOli6?) = -cOldA(cOli1,cOli2,cOli3,cOli4,cOli5,cOli6)
       - 1/20*cOlff(cOli1,cOli2,cOli3,cOli4,cOlk2)*f(cOli5,cOli6,cOlk2)
       - 1/20*cOlff(cOli1,cOli2,cOli3,cOli5,cOlk2)*f(cOli4,cOli6,cOlk2)
       - 1/30*cOlff(cOli1,cOli2,cOli3,cOlk2,cOli4)*f(cOli5,cOli6,cOlk2)
       - 1/30*cOlff(cOli1,cOli2,cOli3,cOlk2,cOli5)*f(cOli4,cOli6,cOlk2)
       - 1/3*cOlff(cOli1,cOli2,cOli3,cOlk2,cOli6)*f(cOli4,cOli5,cOlk2)
       - 1/20*cOlff(cOli1,cOli2,cOli4,cOli3,cOlk2)*f(cOli5,cOli6,cOlk2)
       - 1/20*cOlff(cOli1,cOli2,cOli4,cOli5,cOlk2)*f(cOli3,cOli6,cOlk2)
       - 1/30*cOlff(cOli1,cOli2,cOli4,cOlk2,cOli3)*f(cOli5,cOli6,cOlk2)
       - 1/30*cOlff(cOli1,cOli2,cOli4,cOlk2,cOli5)*f(cOli3,cOli6,cOlk2)
       - 7/30*cOlff(cOli1,cOli2,cOli4,cOlk2,cOli6)*f(cOli3,cOli5,cOlk2)
       - 1/20*cOlff(cOli1,cOli2,cOli5,cOli3,cOlk2)*f(cOli4,cOli6,cOlk2)
       - 1/20*cOlff(cOli1,cOli2,cOli5,cOli4,cOlk2)*f(cOli3,cOli6,cOlk2)
       - 1/30*cOlff(cOli1,cOli2,cOli5,cOlk2,cOli3)*f(cOli4,cOli6,cOlk2)
       - 1/30*cOlff(cOli1,cOli2,cOli5,cOlk2,cOli4)*f(cOli3,cOli6,cOlk2)
       - 1/60*cOlff(cOli1,cOli2,cOlk2,cOli3,cOli4)*f(cOli5,cOli6,cOlk2)
       - 1/60*cOlff(cOli1,cOli2,cOlk2,cOli3,cOli5)*f(cOli4,cOli6,cOlk2)
       - 7/60*cOlff(cOli1,cOli2,cOlk2,cOli3,cOli6)*f(cOli4,cOli5,cOlk2)
       - 1/60*cOlff(cOli1,cOli2,cOlk2,cOli4,cOli3)*f(cOli5,cOli6,cOlk2)
       - 1/60*cOlff(cOli1,cOli2,cOlk2,cOli4,cOli5)*f(cOli3,cOli6,cOlk2)
       - 7/60*cOlff(cOli1,cOli2,cOlk2,cOli4,cOli6)*f(cOli3,cOli5,cOlk2)
       - 1/60*cOlff(cOli1,cOli2,cOlk2,cOli5,cOli3)*f(cOli4,cOli6,cOlk2)
       - 1/60*cOlff(cOli1,cOli2,cOlk2,cOli5,cOli4)*f(cOli3,cOli6,cOlk2)
       - 9/20*cOlff(cOli1,cOli2,cOlk2,cOli5,cOli6)*f(cOli3,cOli4,cOlk2)
       - 1/20*cOlff(cOli1,cOli3,cOli2,cOli4,cOlk2)*f(cOli5,cOli6,cOlk2)
       - 1/20*cOlff(cOli1,cOli3,cOli2,cOli5,cOlk2)*f(cOli4,cOli6,cOlk2)
       - 1/30*cOlff(cOli1,cOli3,cOli2,cOlk2,cOli4)*f(cOli5,cOli6,cOlk2)
       - 1/30*cOlff(cOli1,cOli3,cOli2,cOlk2,cOli5)*f(cOli4,cOli6,cOlk2)
       - 1/20*cOlff(cOli1,cOli3,cOli4,cOli2,cOlk2)*f(cOli5,cOli6,cOlk2)
       - 1/30*cOlff(cOli1,cOli3,cOli4,cOlk2,cOli5)*f(cOli2,cOli6,cOlk2)
       - 1/60*cOlff(cOli1,cOli3,cOli4,cOlk2,cOli6)*f(cOli2,cOli5,cOlk2)
       - 1/20*cOlff(cOli1,cOli3,cOli5,cOli2,cOlk2)*f(cOli4,cOli6,cOlk2)
       - 1/30*cOlff(cOli1,cOli3,cOli5,cOlk2,cOli4)*f(cOli2,cOli6,cOlk2)
       - 1/60*cOlff(cOli1,cOli3,cOli5,cOlk2,cOli6)*f(cOli2,cOli4,cOlk2)
       - 1/60*cOlff(cOli1,cOli3,cOlk2,cOli2,cOli4)*f(cOli5,cOli6,cOlk2)
       - 1/60*cOlff(cOli1,cOli3,cOlk2,cOli2,cOli5)*f(cOli4,cOli6,cOlk2)
       - 1/60*cOlff(cOli1,cOli3,cOlk2,cOli4,cOli5)*f(cOli2,cOli6,cOlk2)
       - 1/12*cOlff(cOli1,cOli3,cOlk2,cOli4,cOli6)*f(cOli2,cOli5,cOlk2)
       - 1/60*cOlff(cOli1,cOli3,cOlk2,cOli5,cOli4)*f(cOli2,cOli6,cOlk2)
       - 1/12*cOlff(cOli1,cOli3,cOlk2,cOli5,cOli6)*f(cOli2,cOli4,cOlk2)
       - 1/20*cOlff(cOli1,cOli4,cOli2,cOli3,cOlk2)*f(cOli5,cOli6,cOlk2)
       - 1/30*cOlff(cOli1,cOli4,cOli2,cOlk2,cOli5)*f(cOli3,cOli6,cOlk2)
       - 1/20*cOlff(cOli1,cOli4,cOli3,cOli2,cOlk2)*f(cOli5,cOli6,cOlk2)
       - 1/30*cOlff(cOli1,cOli4,cOli3,cOlk2,cOli5)*f(cOli2,cOli6,cOlk2)
       - 1/60*cOlff(cOli1,cOli4,cOli3,cOlk2,cOli6)*f(cOli2,cOli5,cOlk2)
       - 1/60*cOlff(cOli1,cOli4,cOli5,cOlk2,cOli6)*f(cOli2,cOli3,cOlk2)
       - 1/60*cOlff(cOli1,cOli4,cOlk2,cOli2,cOli5)*f(cOli3,cOli6,cOlk2)
       - 1/60*cOlff(cOli1,cOli4,cOlk2,cOli3,cOli5)*f(cOli2,cOli6,cOlk2)
       - 1/30*cOlff(cOli1,cOli4,cOlk2,cOli3,cOli6)*f(cOli2,cOli5,cOlk2)
       - 1/12*cOlff(cOli1,cOli4,cOlk2,cOli5,cOli6)*f(cOli2,cOli3,cOlk2)
       - 1/60*cOlff(cOli1,cOli5,cOli3,cOlk2,cOli6)*f(cOli2,cOli4,cOlk2)
       - 1/60*cOlff(cOli1,cOli5,cOli4,cOlk2,cOli6)*f(cOli2,cOli3,cOlk2)
       - 1/30*cOlff(cOli1,cOli5,cOlk2,cOli3,cOli6)*f(cOli2,cOli4,cOlk2)
       - 1/30*cOlff(cOli1,cOli5,cOlk2,cOli4,cOli6)*f(cOli2,cOli3,cOlk2)
       - 1/20*cOlff(cOli1,cOlk2,cOli3,cOli4,cOli6)*f(cOli2,cOli5,cOlk2)
       - 3/20*cOlff(cOli1,cOlk2,cOli3,cOli5,cOli6)*f(cOli2,cOli4,cOlk2)
       - 1/20*cOlff(cOli1,cOlk2,cOli4,cOli3,cOli6)*f(cOli2,cOli5,cOlk2)
       - 3/20*cOlff(cOli1,cOlk2,cOli4,cOli5,cOli6)*f(cOli2,cOli3,cOlk2)
       - 1/20*cOlff(cOli1,cOlk2,cOli5,cOli3,cOli6)*f(cOli2,cOli4,cOlk2)
       - 3/20*cOlff(cOli1,cOlk2,cOli5,cOli4,cOli6)*f(cOli2,cOli3,cOlk2);
#call SORT(adjoint-2-b)
id,once,cOlff(cOli1?,cOli2?,cOli3?,cOli4?,cOli5?) = 
		+1/2*cOldA(cOli2,cOli3,cOli4,cOlk3)*f(cOlk3,cOli1,cOli5)
		-1/2*cOldA(cOli1,cOli4,cOli5,cOlk3)*f(cOlk3,cOli2,cOli3)
		-1/2*cOldA(cOli1,cOli3,cOli5,cOlk3)*f(cOlk3,cOli2,cOli4)
		-1/2*cOldA(cOli1,cOli2,cOli5,cOlk3)*f(cOlk3,cOli3,cOli4)
       + 1/12*f(cOli1,cOlk3,cOlk4)*f(cOli2,cOli3,cOlk3)*f(cOli4,cOli5,cOlk4)*cA
       + 1/12*f(cOli1,cOlk3,cOlk4)*f(cOli2,cOli5,cOlk4)*f(cOli3,cOli4,cOlk3)*cA
       + 1/12*f(cOli1,cOli3,cOlk3)*f(cOli2,cOli4,cOlk4)*f(cOli5,cOlk3,cOlk4)*cA
       - 1/6*f(cOli1,cOli5,cOlk3)*f(cOli2,cOli3,cOlk4)*f(cOli4,cOlk3,cOlk4)*cA
       + 1/12*f(cOli1,cOli5,cOlk3)*f(cOli2,cOli4,cOlk4)*f(cOli3,cOlk3,cOlk4)*cA;
id,once,cOlff(cOli1?,cOli2?,cOli3?,cOli4?) = cOldA(cOli1,cOli2,cOli3,cOli4)
		+cA*f(cOli1,cOli2,cOlk5)*f(cOli4,cOli3,cOlk5)/6
		+cA*f(cOli1,cOlk5,cOli4)*f(cOli3,cOli2,cOlk5)/6;
sum cOlk1,cOlk2,cOlk3,cOlk4,cOlk5;
if ( count(cOlff,1) );
	id  cOlff(?a) = T(?a);
	id  T(cOli1?,cOli2?,?a) = -cOlE0(cOli1,cOli2)*T(?a)-cOlE0(cOlk)*i_*f(cOli1,cOli2,cOlk)*T(?a)/2;
	sum cOlk;
	repeat;
	id  cOlE0(?a)*T(cOli1?,?b) = i_*cOlE0(?a,cOli1)*T(?b)
        +i_*sum_(cOlx,1,nargs_(?a),Bernoulli_(cOlx)*cOlE(cOlx,?a,cOli1))*T(?b);
	repeat;
    	id  cOlE(cOlx?pos_,?a,cOli1?) = distrib_(1,1,cOlEa,cOlEb,?a)*cOlEc(cOlx-1,cOli1,cOlk);
	    id  cOlEa(cOli1?)*cOlEb(?a)*cOlEc(cOlx?,cOli2?,cOli3?) = cOlE(cOlx,?a,cOli3)*i_*f(cOli1,cOli2,cOli3);
    	sum cOlk;
	endrepeat;
	id  cOlE(0,?a) = cOlE0(?a);
	id  T = 1;
	if ( match(T(cOli1?)) );
    	id  cOlE0(?a)*T(cOli1?) = i_*cOldA(?a,cOli1)*cOlacc(1-sign_(nargs_(?a)))/2;
		id	cOlacc(cOlx?) = cOlx;
	elseif ( match(T(cOli1?,cOli2?)) );
    	id  cOlE0(cOli1?)*T(cOli2?,cOli3?) = -i_*cA/2*f(cOli1,cOli2,cOli3);
	    id  cOlE0(cOli1?,cOli2?)*T(cOli3?,cOli4?) = -cOldA(cOli1,cOli2,cOli3,cOli4)
    	    -cA*(f(cOli1,cOli3,cOlk)*f(cOli2,cOli4,cOlk)+f(cOli1,cOli4,cOlk)*f(cOli2,cOli3,cOlk))/12;
	    sum cOlk;
	endif;
	endrepeat;
	id  cOldA(cOli1?) = 0;
	id  cOldA(cOli1?,cOli2?) = cA*d_(cOli1,cOli2);
endif;
#call SORT(adjoint-3)
if ( count(cOldA,1) > 0 );
#do isbb = 1,`RANK'
	if ( count(cOlpA`isbb',1) == 0 );
		id,once,cOldA(?a) = cOldA`isbb'(?a);
		ToVector,cOldA`isbb',cOlpA`isbb';
	endif;
#enddo
endif;
repeat;
	if ( count(cOlff,1) == 0 );
		ReplaceLoop,f,a=3,l=all,outfun=cOlff;
		id	cOlff(cOli1?,cOli2?) = -cA*d_(cOli1,cOli2);
		id	cOlff(cOli1?,cOli2?,cOli3?) = cA/2*f(cOli1,cOli2,cOli3);
	endif;
endrepeat;
if ( count(cOlff,1) ) redefine iStageB "0";
#call SORT(adjoint-4)
#enddo
#endif
#endprocedure
*
*--#] adjoint : 

AutoDeclare Index i,j,k;
CFunction acc;
Symbol x;
.global
G	QQ`SIZE' = <T(i1,i2,j1)>*...*<T(i`SIZE',i{`SIZE'+1},j`SIZE')>
	*<T(k1,k2,j1)>*...*<T(k`SIZE',k{`SIZE'+1},j`SIZE')>
	*replace_(i{`SIZE'+1},i1,k{`SIZE'+1},k1);

G	GG`SIZE' = <f(i1,i2,j1)>*...*<f(i`SIZE',i{`SIZE'+1},j`SIZE')>
	*<f(k1,k2,j1)>*...*<f(k`SIZE',k{`SIZE'+1},j`SIZE')>
	*replace_(i{`SIZE'+1},i1,k{`SIZE'+1},k1);

Sum i1,...,i{`SIZE'*2},j1,...,j`SIZE',k1,...,k`SIZE';
.sort

#define ik1 "0"
repeat id T(cOli1?,cOli2?,?a)*T(cOli2?,cOli3?,?b) = T(cOli1,cOli3,?a,?b);
id  T(cOli1?,cOli1?,?a) = cOlTr(?a);
*
* We work only with the TR. If there are leftover T we don't
* do anything.
*
if ( count(cOlTr,1) > 0 ) redefine ik1 "1";
#call SORT(color-1)
#call adjoint
#if `ik1' != 0
.sort:ik1;
#endif

.end