https://github.com/vermaseren/form/issues/507

507-1.frm

On Highfirst;

*--#[ jodavies_include :

* Paste denexpand.prc

* Expand a "den" function DEN in the variable EP. This is useful in cases
* where we have isolated the dependence on the EP symbol, and cannot merge
* everything into a common PolyRatFun without exceeding MaxTermSize.

* Den is expanded, and the resulting EP-independent coefficients are merged
* into a PolyRatFun PRF as we go. This should most likely be the same name
* as an existing PolyRatFun of EP-independent variables in the input.

Symbol denexpn,denexpm,denexpo,denexpi,denexpj;
CFunction denexpnum,denexpcollect,denexpprftmp;
.sort

#procedure denexpand(DEN,EP,DEPTH,PRF)

	#message denexpand.prc: expanding `DEN' for small `EP', depth `DEPTH'

*	Make sure the sorting mode is highfirst, to work around Issue 336
*	https://github.com/vermaseren/form/issues/336
	On highfirst;

*	PolyRatFun should be disabled initially
	PolyRatFun;

*	Check for multi-argument DEN, this is not expected input
	If (Match(`DEN'(denexpn?,denexpm?)) > 0);
		Print "denexpand.prc: Error: unexpected multi-argument `DEN': %t";
		Exit;
	EndIf;

*	Check PRF doesn't contain the expansion symbol:
	Argument `PRF';
		If (Count(`EP',1) != 0);
			Print "denexpand.prc: Error: unexpected `EP' in `PRF': %t";
			Exit;
		EndIf;
	EndArgument;

*	Extract any poles in EP:
	FactArg `DEN';
	ChainOut `DEN';
	Identify `DEN'(`EP') = 1/`EP';

*	Discard terms which are already beyond DEPTH, now that we are sure DEN contains no poles in EP:
	If (Count(`EP',1) > `DEPTH') Discard;
*	Discard EP powers within DEN, if they will produce terms beyond DEPTH.
*	Account for external negative powers:
	$denexptermEPmax = - count_(`EP',1) + `DEPTH';
	Argument `DEN';
		If (Count(`EP',1) > $denexptermEPmax) Discard;
	EndArgument;

*	Bring DENs to a unique form:
	FactArg `DEN';
	ChainOut `DEN';
	Normalize,^-1 `DEN';
	Identify `DEN'(1) = 1;

*	Record the highest order EP pole
	#$denexpMinEPpole = 0;                                                  
	If (Count(`EP',1) < $denexpMinEPpole) $denexpMinEPpole = count_(`EP',1);

*	Collect together the relevant terms
	AntiBracket `DEN',`EP',`PRF';
	ModuleOption,local $denexptermEPmax;
	ModuleOption,minimum $denexpMinEPpole;
	.sort:DE prep;
*	WORKAROUND (Issue 332)
	#if {-`$denexpMinEPpole'+`DEPTH'} <= 0
		#message denexpand.prc: setting power limit denexpEP(:1);
		Symbol denexpEP(:1);
	#else
		#message denexpand.prc: setting power limit denexpEP(:{-`$denexpMinEPpole'+`DEPTH'});
		Symbol denexpEP(:{-`$denexpMinEPpole'+`DEPTH'});
	#endif


	PolyRatFun `PRF';
	Collect denexpcollect;
*	This sort is necessary to work around Issue 336 ?
	.sort:DE coll;

	#ifdef `DENEXPANDVERBOSE'
		#$term = 1;
		Print "Starting Term %$ in thread %w ...", $term;
		$term = $term + 1;
	#endif

*	Fully compute each coefficient inside a Term environment. Since we used AntiBracket
*	above, the resulting series expansion will never merge with the other terms.
	Term;
*		Open the Collect, so that the PolyRatFuns will properly merge. They do not do so,
*		if we were to attempt the expansion inside the argument of denexpcollect.
		Identify denexpcollect(denexpn?) = denexpn;

*		Convert to power-limited EP symbol for speed, but only external powers and inside
*		the function we are expanding. The symbol might appear in tags, etc.
		Identify `EP' = denexpEP;
		Identify 1/`EP' = 1/denexpEP;
		Argument `DEN';
			Multiply replace_(`EP', denexpEP);
		EndArgument;

*		Change DEN notation
		SplitArg ((denexpEP)) `DEN';
		Transform `DEN' addargs(2,last);
*		A single-argument DEN has no EP dependence and can already be merged with PRF. Typically
*		this doesn't happen, but it depends how the input was prepared.
		Identify `DEN'(denexpj?) = `PRF'(1,denexpj);


*		Loop through the powers of denexpEP which should appear in the output, and fully compute each in turn.
		#do denexpPow = `$denexpMinEPpole',`DEPTH'

			Repeat;
				If (Count(denexpEP,1) == `denexpPow');
*					Expand one DEN by one power, then simplify
					Identify,Once `DEN'(denexpn?,denexpm?) =
						`PRF'(1,denexpn) - `PRF'(1,denexpn)*denexpnum(denexpm/denexpEP)*denexpEP*`DEN'(denexpn,denexpm);

*					Now extract possible denexpEP dependence, without fully expanding denexpnum:
					Repeat;
						SplitArg ((denexpEP)) denexpnum;
						Identify denexpnum(denexpEP) = denexpEP;
						Identify denexpnum(denexpn?) = `PRF'(denexpn,1);
						Transform denexpnum addargs(2,last);
						Identify denexpnum(denexpn?,denexpm?) = `PRF'(denexpn,1) + denexpEP*denexpnum(denexpm/denexpEP);
					EndRepeat;
				EndIf;
			EndRepeat;

*			Before starting the next power, remove terms in the DEN which would produce denexpEP powers beyond DEPTH
			$denexptermEPmax = - count_(denexpEP,1) + `DEPTH'; 
			Argument `DEN';
				If (Count(denexpEP,1) > $denexptermEPmax) Discard;
			EndArgument;
			Identify `DEN'(denexpn?,0) = `PRF'(1,denexpn);

			Sort;
		#enddo

	EndTerm;


*	Replace original EP symbol
	Multiply replace_(denexpEP,`EP');
	#ifdef `DENEXPANDVERBOSE'
		ModuleOption,local $term;
	#endif
	ModuleOption,local $denexptermEPmax;
	.sort:DE exp;
	#printtimes
*	End with PolyRatFun disabled again
	PolyRatFun;

#endprocedure


* Paste coefftools.h

Symbol coefftx,coeffty;

* Convert assumed previous PolyRatFun p (PolyRatFun should now be disabled)
* into products of num and den factors n and d.

#procedure prf2numden(prf,num,den)
	#message coefftools.h: prf2numden(`prf',`num',`den')
	Identify `prf'(coefftx?,coeffty?) = `num'(coefftx)*`den'(coeffty);
	FactArg `num';
	FactArg `den';
	ChainOut `num';
	ChainOut `den';
	Identify `num'(coefftx?number_) = coefftx;
	Identify `den'(coeffty?number_) = 1/coeffty;
#endprocedure


* Convert products of num and den factors n and d into a PolyRatFun p.
* If there are many such factors it is usually faster to loop over moving
* them into the PolyRatFun one at a time.

#procedure numden2prf(num,den,prf)
	#message coefftools.h: numden2prf(`num',`den',`prf')
*	Check for functions inside num, den, which will cause a crash:
	#$termprint = 0;
	Argument `num';
		If (Match(`num'?(?a)));
			$termprint = 1;
		EndIf;
	EndArgument;
	Argument `den';
		If (Match(`num'?(?a)));
			$termprint = 1;
		EndIf;
	EndArgument;
	If ($termprint == 1);
		Print "numden2prf: found a function inside `num' or `den': %t";
		Exit;
	EndIf;
	#do coeffti = 1,1
		#$coeffti = 1;
*		#message coefftools.h: numden2prf loop...
*		Identify,once `num'(coefftx?) = `prf'(coefftx,1);
*		Identify,once `den'(coeffty?) = `prf'(1,coeffty);
		Identify `num'(coefftx?) = `prf'(coefftx,1);
		Identify `den'(coeffty?) = `prf'(1,coeffty);
		If (Count(`num',1,`den',1) > 0) $coeffti = 0;
		ModuleOption,minimum $coeffti;
		ModuleOption,local $termprint;
		.sort(PolyRatFun=`prf')
		#redefine coeffti "`$coeffti'"
		#printtimes
	#enddo
	PolyRatFun;
*	Look for bug Issue 336. It is recommended to always use highfirst sort mode to avoid this.
	If (Match(`prf'(coefftx?,0)) > 0);
		Print "coefftools.h: numden2prf(`num',`den',`prf'): error: division by zero in: %t";
		Exit;
	EndIf;
	.sort
#endprocedure


* Extract S dependence from a numerator function n, which may contain a large sum of terms.
* For eg, #call extractfromnum(ep) will convert num(1+ep*x+ep*y) into 1 + ep*num(x+y).

#procedure extractfromnum(num,S)
	#message coefftools.h: extractfromnum(`num',`S')
	Repeat;
		SplitArg ((`S')) `num';
*		issue-269 workaround, see also https://arxiv.org/abs/math-ph/0010025
		Identify `num'(`S') = `S';
		Transform `num' addargs(2,last);
		Identify `num'(coefftx?,coeffty?) = `num'(coefftx) + `S'*`num'(coeffty/`S');
		FactArg `num';
		ChainOut `num';
		Identify `num'(coefftx?number_) = coefftx;
*		For cases which contained nested num, which already factorized the coefficients of S:
		Repeat Identify `num'(`num'(coefftx?)) = `num'(coefftx);
*		The below leads to horrible performance with lots of expansion followed by re-factorization
*		in the code above. I don't remember why I introduced these lines in the first place, instead
*		to the line above.
*		Argument `num';
*			If (Count(`S',1) == 0) Repeat Identify `num'(coefftx?) = coefftx;
*		EndArgument;
	EndRepeat;
#endprocedure


* It is sometimes useful to rename denominator functions den as deniso, if they
* contain dependence on a particular symbol S.
* For eg, #call isolateden(den,ep,denep) will convert den(x+y)*den(ep*x+y) into
* den(x+y)*denep(ep*x+y).

#procedure isolateden(den,S,deniso)
	#message coefftools.h: isolateden(`den',`S',`deniso')
	FactArg `den';
	ChainOut `den';
	Identify `den'(coefftx?number_) = 1/coefftx;
	SplitArg ((`S')) `den';
	Transform `den' addargs(2,last);
	Identify `den'(`S') = 1/`S';
	Identify `den'(coefftx?,coeffty?) = `deniso'(coefftx+coeffty);
#endprocedure


* Try to factorize the specified list of variables, putting them into products of
* numerator function num. ?S may contain num itself, but should not contain den.
* Try to cancel against any existing den functions.

#procedure mergeterms(num,den,?S)
	AntiBracket `?S';
	.sort
	#message coefftools.h: mergeterms(`num',`den',`?S')
	Collect `num';
	Argument `num';
		Repeat Identify `num'(coefftx?) = coefftx;
		If (Count(`den',1) > 0);
			Print "coefftools.h: mergeterms: `num' should not contain `den' in input: %t";
		EndIf;
	EndArgument;
	FactArg `num';
	ChainOut `num';
	Identify `num'(coefftx?number_) = coefftx;
	Normalize `num';
	Normalize^-1 `den';
	Identify `num'(coefftx?)*`den'(coefftx?) = 1;
#endprocedure
*--#] jodavies_include :

#define NUMTESTS "28"

AutoDeclare Symbol n,m;
Symbol ep;
CFunction prf,den,denep,num,prfexpden;


#do EXPDEPTH = {-1,0,1,5}
	Drop;
	.sort

	Local test1 = prf((n0+ep*n1+ep^2*n2+ep^3*n3),(m0+ep*m1+ep^2*m2+ep^3*m3));
	Local test2 = prf((n0+ep*n1+ep^2*n2+ep^3*n3),ep*(m0+ep*m1+ep^2*m2+ep^3*m3));
	Local test3 = prf(ep*(n0+ep*n1+ep^2*n2+ep^3*n3),(m0+ep*m1+ep^2*m2+ep^3*m3));
	Local test4 = prf(ep*(n0+ep*n1+ep^2*n2+ep^3*n3),ep*(m0+ep*m1+ep^2*m2+ep^3*m3));
	Local test5 = prf((n0+ep*n1+ep^2*n2+ep^3*n3+ep^4*n4+ep^5*n5+ep^6*n6),(m0+ep*m1+ep^2*m2+ep^3*m3));
	Local test6 = prf((n0+ep*n1+ep^2*n2+ep^3*n3+ep^4*n4+ep^5*n5+ep^6*n6),ep*(m0+ep*m1+ep^2*m2+ep^3*m3));
	Local test7 = prf(ep*(n0+ep*n1+ep^2*n2+ep^3*n3+ep^4*n4+ep^5*n5+ep^6*n6),(m0+ep*m1+ep^2*m2+ep^3*m3));
	Local test8 = prf(ep*(n0+ep*n1+ep^2*n2+ep^3*n3+ep^4*n4+ep^5*n5+ep^6*n6),ep*(m0+ep*m1+ep^2*m2+ep^3*m3));
	Local test9 = prf((n0+ep*n1+ep^2*n2+ep^3*n3+ep^4*n4+ep^5*n5+ep^6*n6),(m0+ep*m1+ep^2*m2+ep^3*m3+ep^4*m4+ep^5*m5+ep^6*m6));
	Local test10 = prf((n0+ep*n1+ep^2*n2+ep^3*n3+ep^4*n4+ep^5*n5+ep^6*n6),ep*(m0+ep*m1+ep^2*m2+ep^3*m3+ep^4*m4+ep^5*m5+ep^6*m6));
	Local test11 = prf(ep*(n0+ep*n1+ep^2*n2+ep^3*n3+ep^4*n4+ep^5*n5+ep^6*n6),(m0+ep*m1+ep^2*m2+ep^3*m3+ep^4*m4+ep^5*m5+ep^6*m6));
	Local test12 = prf(ep*(n0+ep*n1+ep^2*n2+ep^3*n3+ep^4*n4+ep^5*n5+ep^6*n6),ep*(m0+ep*m1+ep^2*m2+ep^3*m3+ep^4*m4+ep^5*m5+ep^6*m6));
	Local test13 = prf((n0+ep*n1+ep^2*n2+ep^3*n3),(ep*m0+ep*m1+ep^2*m2+ep^3*m3));
	Local test14 = prf((n0+ep*n1+ep^2*n2+ep^3*n3),ep*(ep*m0+ep*m1+ep^2*m2+ep^3*m3));
	Local test15 = prf(ep*(n0+ep*n1+ep^2*n2+ep^3*n3),(ep*m0+ep*m1+ep^2*m2+ep^3*m3));
	Local test16 = prf(ep*(n0+ep*n1+ep^2*n2+ep^3*n3),ep*(ep*m0+ep*m1+ep^2*m2+ep^3*m3));
	Local test17 = prf(1,1+ep);
	Local test18 = prf(1,ep^5*(1+ep));
	Local test19 = prf(1,ep^5*(1+2*ep));
	Local test20 = prf(1,1);
	Local test21 = prf(1,1+ep)+prf(-1,1)+prf(ep,1+ep);
	Local test22 = prf(1,ep);
	Local test23 = prf(ep,ep);
	Local test24 = prf(0,ep);
	Local test25 = prf(0,1);
	Local test26 = prf(ep^6,1);
	Local test27 = prf(ep^7,ep);
*	This test breaks if the internal power limit is set to zero:
	Local test28 = prf(1,m0+m1+ep*(m0+m1));

*	Format for denexpand
	#call prf2numden(prf,num,den)
	#call extractfromnum(num,ep)
	#call isolateden(den,ep,denep)
	#call numden2prf(num,den,prf)
	.sort

	#call denexpand(denep,ep,`EXPDEPTH',prf)

	Local result1 = n0*prfexpden(m0) + ep*(-(m1*n0) + m0*n1)*prfexpden(m0)^2 + ep^2*(m1^2*n0 - m0*m2*n0 - m0*m1*n1 + m0^2*n2)*prfexpden(m0)^3 + ep^3*(-(m1^3*n0) + 2*m0*m1*m2*n0 - m0^2*m3*n0 + m0*m1^2*n1 - m0^2*m2*n1 - m0^2*m1*n2 + m0^3*n3)*prfexpden(m0)^4 + ep^4*(m1^4*n0 - 3*m0*m1^2*m2*n0 + m0^2*m2^2*n0 + 2*m0^2*m1*m3*n0 - m0*m1^3*n1 + 2*m0^2*m1*m2*n1 - m0^3*m3*n1 + m0^2*m1^2*n2 - m0^3*m2*n2 - m0^3*m1*n3)*prfexpden(m0)^5 + ep^5*(-(m1^5*n0) + 4*m0*m1^3*m2*n0 - 3*m0^2*m1*m2^2*n0 - 3*m0^2*m1^2*m3*n0 + 2*m0^3*m2*m3*n0 + m0*m1^4*n1 - 3*m0^2*m1^2*m2*n1 + m0^3*m2^2*n1 + 2*m0^3*m1*m3*n1 - m0^2*m1^3*n2 + 2*m0^3*m1*m2*n2 - m0^4*m3*n2 + m0^3*m1^2*n3 - m0^4*m2*n3)*prfexpden(m0)^6;
	Local result2 = n0*prfexpden(ep)*prfexpden(m0) + (-(m1*n0) + m0*n1)*prfexpden(m0)^2 + ep*(m1^2*n0 - m0*m2*n0 - m0*m1*n1 + m0^2*n2)*prfexpden(m0)^3 + ep^2*(n3*prfexpden(m0) - m1*n2*prfexpden(m0)^2 + (m1^2 - m0*m2)*n1*prfexpden(m0)^3 + (-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n0*prfexpden(m0)^4) + ep^3*(-(m1*n3*prfexpden(m0)^2) + (m1^2 - m0*m2)*n2*prfexpden(m0)^3 + (-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n1*prfexpden(m0)^4 + (m1^4 - 3*m0*m1^2*m2 + m0^2*m2^2 + 2*m0^2*m1*m3)*n0*prfexpden(m0)^5) + ep^4*((m1^2 - m0*m2)*n3*prfexpden(m0)^3 + (-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n2*prfexpden(m0)^4 + (m1^4 - 3*m0*m1^2*m2 + m0^2*m2^2 + 2*m0^2*m1*m3)*n1*prfexpden(m0)^5 + (-m1^5 + 4*m0*m1^3*m2 - 3*m0^2*m1*m2^2 - 3*m0^2*m1^2*m3 + 2*m0^3*m2*m3)*n0*prfexpden(m0)^6) + ep^5*((-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n3*prfexpden(m0)^4 + (m1^4 - 3*m0*m1^2*m2 + m0^2*m2^2 + 2*m0^2*m1*m3)*n2*prfexpden(m0)^5 + (-m1^5 + 4*m0*m1^3*m2 - 3*m0^2*m1*m2^2 - 3*m0^2*m1^2*m3 + 2*m0^3*m2*m3)*n1*prfexpden(m0)^6 + (m1^6 - 5*m0*m1^4*m2 + 6*m0^2*m1^2*m2^2 - m0^3*m2^3 + 4*m0^2*m1^3*m3 - 6*m0^3*m1*m2*m3 + m0^4*m3^2)*n0*prfexpden(m0)^7);
	Local result3 = ep*n0*prfexpden(m0) + ep^2*(-(m1*n0) + m0*n1)*prfexpden(m0)^2 + ep^3*(m1^2*n0 - m0*m2*n0 - m0*m1*n1 + m0^2*n2)*prfexpden(m0)^3 + ep^4*(-(m1^3*n0) + 2*m0*m1*m2*n0 - m0^2*m3*n0 + m0*m1^2*n1 - m0^2*m2*n1 - m0^2*m1*n2 + m0^3*n3)*prfexpden(m0)^4 + ep^5*(m1^4*n0 - 3*m0*m1^2*m2*n0 + m0^2*m2^2*n0 + 2*m0^2*m1*m3*n0 - m0*m1^3*n1 + 2*m0^2*m1*m2*n1 - m0^3*m3*n1 + m0^2*m1^2*n2 - m0^3*m2*n2 - m0^3*m1*n3)*prfexpden(m0)^5;
	Local result4 = n0*prfexpden(m0) + ep*(-(m1*n0) + m0*n1)*prfexpden(m0)^2 + ep^2*(m1^2*n0 - m0*m2*n0 - m0*m1*n1 + m0^2*n2)*prfexpden(m0)^3 + ep^3*(-(m1^3*n0) + 2*m0*m1*m2*n0 - m0^2*m3*n0 + m0*m1^2*n1 - m0^2*m2*n1 - m0^2*m1*n2 + m0^3*n3)*prfexpden(m0)^4 + ep^4*(m1^4*n0 - 3*m0*m1^2*m2*n0 + m0^2*m2^2*n0 + 2*m0^2*m1*m3*n0 - m0*m1^3*n1 + 2*m0^2*m1*m2*n1 - m0^3*m3*n1 + m0^2*m1^2*n2 - m0^3*m2*n2 - m0^3*m1*n3)*prfexpden(m0)^5 + ep^5*(-(m1^5*n0) + 4*m0*m1^3*m2*n0 - 3*m0^2*m1*m2^2*n0 - 3*m0^2*m1^2*m3*n0 + 2*m0^3*m2*m3*n0 + m0*m1^4*n1 - 3*m0^2*m1^2*m2*n1 + m0^3*m2^2*n1 + 2*m0^3*m1*m3*n1 - m0^2*m1^3*n2 + 2*m0^3*m1*m2*n2 - m0^4*m3*n2 + m0^3*m1^2*n3 - m0^4*m2*n3)*prfexpden(m0)^6;
	Local result5 = n0*prfexpden(m0) + ep*(-(m1*n0) + m0*n1)*prfexpden(m0)^2 + ep^2*(m1^2*n0 - m0*m2*n0 - m0*m1*n1 + m0^2*n2)*prfexpden(m0)^3 + ep^3*(-(m1^3*n0) + 2*m0*m1*m2*n0 - m0^2*m3*n0 + m0*m1^2*n1 - m0^2*m2*n1 - m0^2*m1*n2 + m0^3*n3)*prfexpden(m0)^4 + ep^4*(m1^4*n0 - 3*m0*m1^2*m2*n0 + m0^2*m2^2*n0 + 2*m0^2*m1*m3*n0 - m0*m1^3*n1 + 2*m0^2*m1*m2*n1 - m0^3*m3*n1 + m0^2*m1^2*n2 - m0^3*m2*n2 - m0^3*m1*n3 + m0^4*n4)*prfexpden(m0)^5 + ep^5*(-(m1^5*n0) + 4*m0*m1^3*m2*n0 - 3*m0^2*m1*m2^2*n0 - 3*m0^2*m1^2*m3*n0 + 2*m0^3*m2*m3*n0 + m0*m1^4*n1 - 3*m0^2*m1^2*m2*n1 + m0^3*m2^2*n1 + 2*m0^3*m1*m3*n1 - m0^2*m1^3*n2 + 2*m0^3*m1*m2*n2 - m0^4*m3*n2 + m0^3*m1^2*n3 - m0^4*m2*n3 - m0^4*m1*n4 + m0^5*n5)*prfexpden(m0)^6;
	Local result6 = n0*prfexpden(ep)*prfexpden(m0) + (-(m1*n0) + m0*n1)*prfexpden(m0)^2 + ep*(m1^2*n0 - m0*m2*n0 - m0*m1*n1 + m0^2*n2)*prfexpden(m0)^3 + ep^2*(n3*prfexpden(m0) - m1*n2*prfexpden(m0)^2 + (m1^2 - m0*m2)*n1*prfexpden(m0)^3 + (-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n0*prfexpden(m0)^4) + ep^3*(n4*prfexpden(m0) - m1*n3*prfexpden(m0)^2 + (m1^2 - m0*m2)*n2*prfexpden(m0)^3 + (-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n1*prfexpden(m0)^4 + (m1^4 - 3*m0*m1^2*m2 + m0^2*m2^2 + 2*m0^2*m1*m3)*n0*prfexpden(m0)^5) + ep^4*(n5*prfexpden(m0) - m1*n4*prfexpden(m0)^2 + (m1^2 - m0*m2)*n3*prfexpden(m0)^3 + (-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n2*prfexpden(m0)^4 + (m1^4 - 3*m0*m1^2*m2 + m0^2*m2^2 + 2*m0^2*m1*m3)*n1*prfexpden(m0)^5 + (-m1^5 + 4*m0*m1^3*m2 - 3*m0^2*m1*m2^2 - 3*m0^2*m1^2*m3 + 2*m0^3*m2*m3)*n0*prfexpden(m0)^6) + ep^5*(n6*prfexpden(m0) - m1*n5*prfexpden(m0)^2 + (m1^2 - m0*m2)*n4*prfexpden(m0)^3 + (-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n3*prfexpden(m0)^4 + (m1^4 - 3*m0*m1^2*m2 + m0^2*m2^2 + 2*m0^2*m1*m3)*n2*prfexpden(m0)^5 + (-m1^5 + 4*m0*m1^3*m2 - 3*m0^2*m1*m2^2 - 3*m0^2*m1^2*m3 + 2*m0^3*m2*m3)*n1*prfexpden(m0)^6 + (m1^6 - 5*m0*m1^4*m2 + 6*m0^2*m1^2*m2^2 - m0^3*m2^3 + 4*m0^2*m1^3*m3 - 6*m0^3*m1*m2*m3 + m0^4*m3^2)*n0*prfexpden(m0)^7);
	Local result7 = ep*n0*prfexpden(m0) + ep^2*(-(m1*n0) + m0*n1)*prfexpden(m0)^2 + ep^3*(m1^2*n0 - m0*m2*n0 - m0*m1*n1 + m0^2*n2)*prfexpden(m0)^3 + ep^4*(-(m1^3*n0) + 2*m0*m1*m2*n0 - m0^2*m3*n0 + m0*m1^2*n1 - m0^2*m2*n1 - m0^2*m1*n2 + m0^3*n3)*prfexpden(m0)^4 + ep^5*(m1^4*n0 - 3*m0*m1^2*m2*n0 + m0^2*m2^2*n0 + 2*m0^2*m1*m3*n0 - m0*m1^3*n1 + 2*m0^2*m1*m2*n1 - m0^3*m3*n1 + m0^2*m1^2*n2 - m0^3*m2*n2 - m0^3*m1*n3 + m0^4*n4)*prfexpden(m0)^5;
	Local result8 = n0*prfexpden(m0) + ep*(-(m1*n0) + m0*n1)*prfexpden(m0)^2 + ep^2*(m1^2*n0 - m0*m2*n0 - m0*m1*n1 + m0^2*n2)*prfexpden(m0)^3 + ep^3*(-(m1^3*n0) + 2*m0*m1*m2*n0 - m0^2*m3*n0 + m0*m1^2*n1 - m0^2*m2*n1 - m0^2*m1*n2 + m0^3*n3)*prfexpden(m0)^4 + ep^4*(m1^4*n0 - 3*m0*m1^2*m2*n0 + m0^2*m2^2*n0 + 2*m0^2*m1*m3*n0 - m0*m1^3*n1 + 2*m0^2*m1*m2*n1 - m0^3*m3*n1 + m0^2*m1^2*n2 - m0^3*m2*n2 - m0^3*m1*n3 + m0^4*n4)*prfexpden(m0)^5 + ep^5*(-(m1^5*n0) + 4*m0*m1^3*m2*n0 - 3*m0^2*m1*m2^2*n0 - 3*m0^2*m1^2*m3*n0 + 2*m0^3*m2*m3*n0 + m0*m1^4*n1 - 3*m0^2*m1^2*m2*n1 + m0^3*m2^2*n1 + 2*m0^3*m1*m3*n1 - m0^2*m1^3*n2 + 2*m0^3*m1*m2*n2 - m0^4*m3*n2 + m0^3*m1^2*n3 - m0^4*m2*n3 - m0^4*m1*n4 + m0^5*n5)*prfexpden(m0)^6;
	Local result9 = n0*prfexpden(m0) + ep*(-(m1*n0) + m0*n1)*prfexpden(m0)^2 + ep^2*(m1^2*n0 - m0*m2*n0 - m0*m1*n1 + m0^2*n2)*prfexpden(m0)^3 + ep^3*(n3*prfexpden(m0) - m1*n2*prfexpden(m0)^2 + (m1^2 - m0*m2)*n1*prfexpden(m0)^3 + (-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n0*prfexpden(m0)^4) + ep^4*(n4*prfexpden(m0) - m1*n3*prfexpden(m0)^2 + (m1^2 - m0*m2)*n2*prfexpden(m0)^3 + (-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n1*prfexpden(m0)^4 + (m1^4 - 3*m0*m1^2*m2 + m0^2*m2^2 + 2*m0^2*m1*m3 - m0^3*m4)*n0*prfexpden(m0)^5) + ep^5*(n5*prfexpden(m0) - m1*n4*prfexpden(m0)^2 + (m1^2 - m0*m2)*n3*prfexpden(m0)^3 + (-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n2*prfexpden(m0)^4 + (m1^4 - 3*m0*m1^2*m2 + m0^2*m2^2 + 2*m0^2*m1*m3 - m0^3*m4)*n1*prfexpden(m0)^5 + (-m1^5 + 4*m0*m1^3*m2 - 3*m0^2*m1*m2^2 - 3*m0^2*m1^2*m3 + 2*m0^3*m2*m3 + 2*m0^3*m1*m4 - m0^4*m5)*n0*prfexpden(m0)^6);
	Local result10 = n0*prfexpden(ep)*prfexpden(m0) + (-(m1*n0) + m0*n1)*prfexpden(m0)^2 + ep*(m1^2*n0 - m0*m2*n0 - m0*m1*n1 + m0^2*n2)*prfexpden(m0)^3 + ep^2*(n3*prfexpden(m0) - m1*n2*prfexpden(m0)^2 + (m1^2 - m0*m2)*n1*prfexpden(m0)^3 + (-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n0*prfexpden(m0)^4) + ep^3*(n4*prfexpden(m0) - m1*n3*prfexpden(m0)^2 + (m1^2 - m0*m2)*n2*prfexpden(m0)^3 + (-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n1*prfexpden(m0)^4 + (m1^4 - 3*m0*m1^2*m2 + m0^2*m2^2 + 2*m0^2*m1*m3 - m0^3*m4)*n0*prfexpden(m0)^5) + ep^4*(n5*prfexpden(m0) - m1*n4*prfexpden(m0)^2 + (m1^2 - m0*m2)*n3*prfexpden(m0)^3 + (-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n2*prfexpden(m0)^4 + (m1^4 - 3*m0*m1^2*m2 + m0^2*m2^2 + 2*m0^2*m1*m3 - m0^3*m4)*n1*prfexpden(m0)^5 + (-m1^5 + 4*m0*m1^3*m2 - 3*m0^2*m1*m2^2 - 3*m0^2*m1^2*m3 + 2*m0^3*m2*m3 + 2*m0^3*m1*m4 - m0^4*m5)*n0*prfexpden(m0)^6) + ep^5*(n6*prfexpden(m0) - m1*n5*prfexpden(m0)^2 + (m1^2 - m0*m2)*n4*prfexpden(m0)^3 + (-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n3*prfexpden(m0)^4 + (m1^4 - 3*m0*m1^2*m2 + m0^2*m2^2 + 2*m0^2*m1*m3 - m0^3*m4)*n2*prfexpden(m0)^5 + (-m1^5 + 4*m0*m1^3*m2 - 3*m0^2*m1*m2^2 - 3*m0^2*m1^2*m3 + 2*m0^3*m2*m3 + 2*m0^3*m1*m4 - m0^4*m5)*n1*prfexpden(m0)^6 + (m1^6 - 5*m0*m1^4*m2 + 6*m0^2*m1^2*m2^2 - m0^3*m2^3 + 4*m0^2*m1^3*m3 - 6*m0^3*m1*m2*m3 + m0^4*m3^2 - 3*m0^3*m1^2*m4 + 2*m0^4*m2*m4 + 2*m0^4*m1*m5 - m0^5*m6)*n0*prfexpden(m0)^7);
	Local result11 = ep*n0*prfexpden(m0) + ep^2*(-(m1*n0) + m0*n1)*prfexpden(m0)^2 + ep^3*(m1^2*n0 - m0*m2*n0 - m0*m1*n1 + m0^2*n2)*prfexpden(m0)^3 + ep^4*(-(m1^3*n0) + 2*m0*m1*m2*n0 - m0^2*m3*n0 + m0*m1^2*n1 - m0^2*m2*n1 - m0^2*m1*n2 + m0^3*n3)*prfexpden(m0)^4 + ep^5*(m1^4*n0 - 3*m0*m1^2*m2*n0 + m0^2*m2^2*n0 + 2*m0^2*m1*m3*n0 - m0^3*m4*n0 - m0*m1^3*n1 + 2*m0^2*m1*m2*n1 - m0^3*m3*n1 + m0^2*m1^2*n2 - m0^3*m2*n2 - m0^3*m1*n3 + m0^4*n4)*prfexpden(m0)^5;
	Local result12 = n0*prfexpden(m0) + ep*(-(m1*n0) + m0*n1)*prfexpden(m0)^2 + ep^2*(m1^2*n0 - m0*m2*n0 - m0*m1*n1 + m0^2*n2)*prfexpden(m0)^3 + ep^3*(n3*prfexpden(m0) - m1*n2*prfexpden(m0)^2 + (m1^2 - m0*m2)*n1*prfexpden(m0)^3 + (-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n0*prfexpden(m0)^4) + ep^4*(n4*prfexpden(m0) - m1*n3*prfexpden(m0)^2 + (m1^2 - m0*m2)*n2*prfexpden(m0)^3 + (-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n1*prfexpden(m0)^4 + (m1^4 - 3*m0*m1^2*m2 + m0^2*m2^2 + 2*m0^2*m1*m3 - m0^3*m4)*n0*prfexpden(m0)^5) + ep^5*(n5*prfexpden(m0) - m1*n4*prfexpden(m0)^2 + (m1^2 - m0*m2)*n3*prfexpden(m0)^3 + (-m1^3 + 2*m0*m1*m2 - m0^2*m3)*n2*prfexpden(m0)^4 + (m1^4 - 3*m0*m1^2*m2 + m0^2*m2^2 + 2*m0^2*m1*m3 - m0^3*m4)*n1*prfexpden(m0)^5 + (-m1^5 + 4*m0*m1^3*m2 - 3*m0^2*m1*m2^2 - 3*m0^2*m1^2*m3 + 2*m0^3*m2*m3 + 2*m0^3*m1*m4 - m0^4*m5)*n0*prfexpden(m0)^6);
	Local result13 = n1*prfexpden(m0 + m1) + n0*prfexpden(ep)*prfexpden(m0 + m1) - m2*n0*prfexpden(m0 + m1)^2 + ep*(n2*prfexpden(m0 + m1) - m2*n1*prfexpden(m0 + m1)^2 + (m2^2 - m0*m3 - m1*m3)*n0*prfexpden(m0 + m1)^3) + ep^2*(n3*prfexpden(m0 + m1) - m2*n2*prfexpden(m0 + m1)^2 + (m2^2 - m0*m3 - m1*m3)*n1*prfexpden(m0 + m1)^3 + (-m2^3 + 2*m0*m2*m3 + 2*m1*m2*m3)*n0*prfexpden(m0 + m1)^4) + ep^3*(-(m2*n3*prfexpden(m0 + m1)^2) + (m2^2 - m0*m3 - m1*m3)*n2*prfexpden(m0 + m1)^3 + (-m2^3 + 2*m0*m2*m3 + 2*m1*m2*m3)*n1*prfexpden(m0 + m1)^4 + n0*prfexpden(m0 + m1)*(m3^2*prfexpden(m0 + m1)^2 - 3*m2^2*m3*prfexpden(m0 + m1)^3 + m2^4*prfexpden(m0 + m1)^4)) + ep^4*((m2^2 - m0*m3 - m1*m3)*n3*prfexpden(m0 + m1)^3 + (-m2^3 + 2*m0*m2*m3 + 2*m1*m2*m3)*n2*prfexpden(m0 + m1)^4 + n1*prfexpden(m0 + m1)*(m3^2*prfexpden(m0 + m1)^2 - 3*m2^2*m3*prfexpden(m0 + m1)^3 + m2^4*prfexpden(m0 + m1)^4) + n0*prfexpden(m0 + m1)*(-3*m2*m3^2*prfexpden(m0 + m1)^3 + 4*m2^3*m3*prfexpden(m0 + m1)^4 - m2^5*prfexpden(m0 + m1)^5)) + ep^5*((-m2^3 + 2*m0*m2*m3 + 2*m1*m2*m3)*n3*prfexpden(m0 + m1)^4 + n2*prfexpden(m0 + m1)*(m3^2*prfexpden(m0 + m1)^2 - 3*m2^2*m3*prfexpden(m0 + m1)^3 + m2^4*prfexpden(m0 + m1)^4) + n1*prfexpden(m0 + m1)*(-3*m2*m3^2*prfexpden(m0 + m1)^3 + 4*m2^3*m3*prfexpden(m0 + m1)^4 - m2^5*prfexpden(m0 + m1)^5) + n0*prfexpden(m0 + m1)*(-(m3^3*prfexpden(m0 + m1)^3) + 6*m2^2*m3^2*prfexpden(m0 + m1)^4 - 5*m2^4*m3*prfexpden(m0 + m1)^5 + m2^6*prfexpden(m0 + m1)^6));
	Local result14 = n2*prfexpden(m0 + m1) + n0*prfexpden(ep)^2*prfexpden(m0 + m1) - m2*n1*prfexpden(m0 + m1)^2 + (m2^2 - m0*m3 - m1*m3)*n0*prfexpden(m0 + m1)^3 + prfexpden(ep)*(n1*prfexpden(m0 + m1) - m2*n0*prfexpden(m0 + m1)^2) + ep*(n3*prfexpden(m0 + m1) - m2*n2*prfexpden(m0 + m1)^2 + (m2^2 - m0*m3 - m1*m3)*n1*prfexpden(m0 + m1)^3 + (-m2^3 + 2*m0*m2*m3 + 2*m1*m2*m3)*n0*prfexpden(m0 + m1)^4) + ep^2*(-(m2*n3*prfexpden(m0 + m1)^2) + (m2^2 - m0*m3 - m1*m3)*n2*prfexpden(m0 + m1)^3 + (-m2^3 + 2*m0*m2*m3 + 2*m1*m2*m3)*n1*prfexpden(m0 + m1)^4 + n0*prfexpden(m0 + m1)*(m3^2*prfexpden(m0 + m1)^2 - 3*m2^2*m3*prfexpden(m0 + m1)^3 + m2^4*prfexpden(m0 + m1)^4)) + ep^3*((m2^2 - m0*m3 - m1*m3)*n3*prfexpden(m0 + m1)^3 + (-m2^3 + 2*m0*m2*m3 + 2*m1*m2*m3)*n2*prfexpden(m0 + m1)^4 + n1*prfexpden(m0 + m1)*(m3^2*prfexpden(m0 + m1)^2 - 3*m2^2*m3*prfexpden(m0 + m1)^3 + m2^4*prfexpden(m0 + m1)^4) + n0*prfexpden(m0 + m1)*(-3*m2*m3^2*prfexpden(m0 + m1)^3 + 4*m2^3*m3*prfexpden(m0 + m1)^4 - m2^5*prfexpden(m0 + m1)^5)) + ep^4*((-m2^3 + 2*m0*m2*m3 + 2*m1*m2*m3)*n3*prfexpden(m0 + m1)^4 + n2*prfexpden(m0 + m1)*(m3^2*prfexpden(m0 + m1)^2 - 3*m2^2*m3*prfexpden(m0 + m1)^3 + m2^4*prfexpden(m0 + m1)^4) + n1*prfexpden(m0 + m1)*(-3*m2*m3^2*prfexpden(m0 + m1)^3 + 4*m2^3*m3*prfexpden(m0 + m1)^4 - m2^5*prfexpden(m0 + m1)^5) + n0*prfexpden(m0 + m1)*(-(m3^3*prfexpden(m0 + m1)^3) + 6*m2^2*m3^2*prfexpden(m0 + m1)^4 - 5*m2^4*m3*prfexpden(m0 + m1)^5 + m2^6*prfexpden(m0 + m1)^6)) + ep^5*(n3*prfexpden(m0 + m1)*(m3^2*prfexpden(m0 + m1)^2 - 3*m2^2*m3*prfexpden(m0 + m1)^3 + m2^4*prfexpden(m0 + m1)^4) + n2*prfexpden(m0 + m1)*(-3*m2*m3^2*prfexpden(m0 + m1)^3 + 4*m2^3*m3*prfexpden(m0 + m1)^4 - m2^5*prfexpden(m0 + m1)^5) + n1*prfexpden(m0 + m1)*(-(m3^3*prfexpden(m0 + m1)^3) + 6*m2^2*m3^2*prfexpden(m0 + m1)^4 - 5*m2^4*m3*prfexpden(m0 + m1)^5 + m2^6*prfexpden(m0 + m1)^6) + n0*prfexpden(m0 + m1)*(4*m2*m3^3*prfexpden(m0 + m1)^4 - 10*m2^3*m3^2*prfexpden(m0 + m1)^5 + 6*m2^5*m3*prfexpden(m0 + m1)^6 - m2^7*prfexpden(m0 + m1)^7));
	Local result15 = n0*prfexpden(m0 + m1) + ep*(n1*prfexpden(m0 + m1) - m2*n0*prfexpden(m0 + m1)^2) + ep^2*(n2*prfexpden(m0 + m1) - m2*n1*prfexpden(m0 + m1)^2 + (m2^2 - m0*m3 - m1*m3)*n0*prfexpden(m0 + m1)^3) + ep^3*(n3*prfexpden(m0 + m1) - m2*n2*prfexpden(m0 + m1)^2 + (m2^2 - m0*m3 - m1*m3)*n1*prfexpden(m0 + m1)^3 + (-m2^3 + 2*m0*m2*m3 + 2*m1*m2*m3)*n0*prfexpden(m0 + m1)^4) + ep^4*(-(m2*n3*prfexpden(m0 + m1)^2) + (m2^2 - m0*m3 - m1*m3)*n2*prfexpden(m0 + m1)^3 + (-m2^3 + 2*m0*m2*m3 + 2*m1*m2*m3)*n1*prfexpden(m0 + m1)^4 + n0*prfexpden(m0 + m1)*(m3^2*prfexpden(m0 + m1)^2 - 3*m2^2*m3*prfexpden(m0 + m1)^3 + m2^4*prfexpden(m0 + m1)^4)) + ep^5*((m2^2 - m0*m3 - m1*m3)*n3*prfexpden(m0 + m1)^3 + (-m2^3 + 2*m0*m2*m3 + 2*m1*m2*m3)*n2*prfexpden(m0 + m1)^4 + n1*prfexpden(m0 + m1)*(m3^2*prfexpden(m0 + m1)^2 - 3*m2^2*m3*prfexpden(m0 + m1)^3 + m2^4*prfexpden(m0 + m1)^4) + n0*prfexpden(m0 + m1)*(-3*m2*m3^2*prfexpden(m0 + m1)^3 + 4*m2^3*m3*prfexpden(m0 + m1)^4 - m2^5*prfexpden(m0 + m1)^5));
	Local result16 = n1*prfexpden(m0 + m1) + n0*prfexpden(ep)*prfexpden(m0 + m1) - m2*n0*prfexpden(m0 + m1)^2 + ep*(n2*prfexpden(m0 + m1) - m2*n1*prfexpden(m0 + m1)^2 + (m2^2 - m0*m3 - m1*m3)*n0*prfexpden(m0 + m1)^3) + ep^2*(n3*prfexpden(m0 + m1) - m2*n2*prfexpden(m0 + m1)^2 + (m2^2 - m0*m3 - m1*m3)*n1*prfexpden(m0 + m1)^3 + (-m2^3 + 2*m0*m2*m3 + 2*m1*m2*m3)*n0*prfexpden(m0 + m1)^4) + ep^3*(-(m2*n3*prfexpden(m0 + m1)^2) + (m2^2 - m0*m3 - m1*m3)*n2*prfexpden(m0 + m1)^3 + (-m2^3 + 2*m0*m2*m3 + 2*m1*m2*m3)*n1*prfexpden(m0 + m1)^4 + n0*prfexpden(m0 + m1)*(m3^2*prfexpden(m0 + m1)^2 - 3*m2^2*m3*prfexpden(m0 + m1)^3 + m2^4*prfexpden(m0 + m1)^4)) + ep^4*((m2^2 - m0*m3 - m1*m3)*n3*prfexpden(m0 + m1)^3 + (-m2^3 + 2*m0*m2*m3 + 2*m1*m2*m3)*n2*prfexpden(m0 + m1)^4 + n1*prfexpden(m0 + m1)*(m3^2*prfexpden(m0 + m1)^2 - 3*m2^2*m3*prfexpden(m0 + m1)^3 + m2^4*prfexpden(m0 + m1)^4) + n0*prfexpden(m0 + m1)*(-3*m2*m3^2*prfexpden(m0 + m1)^3 + 4*m2^3*m3*prfexpden(m0 + m1)^4 - m2^5*prfexpden(m0 + m1)^5)) + ep^5*((-m2^3 + 2*m0*m2*m3 + 2*m1*m2*m3)*n3*prfexpden(m0 + m1)^4 + n2*prfexpden(m0 + m1)*(m3^2*prfexpden(m0 + m1)^2 - 3*m2^2*m3*prfexpden(m0 + m1)^3 + m2^4*prfexpden(m0 + m1)^4) + n1*prfexpden(m0 + m1)*(-3*m2*m3^2*prfexpden(m0 + m1)^3 + 4*m2^3*m3*prfexpden(m0 + m1)^4 - m2^5*prfexpden(m0 + m1)^5) + n0*prfexpden(m0 + m1)*(-(m3^3*prfexpden(m0 + m1)^3) + 6*m2^2*m3^2*prfexpden(m0 + m1)^4 - 5*m2^4*m3*prfexpden(m0 + m1)^5 + m2^6*prfexpden(m0 + m1)^6));
	Local result17 = 1 - ep + ep^2 - ep^3 + ep^4 - ep^5;
	Local result18 = -1 + ep - ep^2 + ep^3 - ep^4 + ep^5 + prfexpden(ep) - prfexpden(ep)^2 + prfexpden(ep)^3 - prfexpden(ep)^4 + prfexpden(ep)^5;
	Local result19 = -32 + 64*ep - 128*ep^2 + 256*ep^3 - 512*ep^4 + 1024*ep^5 + 16*prfexpden(ep) - 8*prfexpden(ep)^2 + 4*prfexpden(ep)^3 - 2*prfexpden(ep)^4 + prfexpden(ep)^5;
	Local result20 = 1;
	Local result21 = 0;
	Local result22 = ep^-1;
	Local result23 = 1;
	Local result24 = 0;
	Local result25 = 0;
	Local result26 = 0;
	Local result27 = 0;
	Local result28 = prfexpden(m0 + m1) - ep*prfexpden(m0 + m1) + ep^2*prfexpden(m0 + m1) - ep^3*prfexpden(m0 + m1) + ep^4*prfexpden(m0 + m1) - ep^5*prfexpden(m0 + m1);
	.sort
	InExpression result1,...,result`NUMTESTS';
		Identify prfexpden(ep) = 1/ep;
		If (Count(ep,1) > `EXPDEPTH') Discard;
	EndInExpression;
	.sort

	PolyRatFun prf;
	Drop;
*	Compare
	#do i = 1,`NUMTESTS'
		Local diff`i' = test`i' - result`i';
	#enddo
*	Normalize notation between FORM and Mathematica results
	Identify prfexpden(n?) = prf(1,n);
	#do i = 0,10
		Identify n`i' = prf(n`i',1);
		Identify m`i' = prf(m`i',1);
	#enddo
	.sort

	#do i = 1,`NUMTESTS'
		#if `ZERO_diff`i'' == 0
			#message Error in test`i' EXPDEPTH `EXPDEPTH'
			#terminate
		#endif
	#enddo

#enddo

.end

507-2.frm

AutoDeclare Symbol n,m;
Symbol ep;
CFunction den,num;
Symbol coefftx,coeffty;

Local test4 = + den(ep)*den(ep^3*m3 + ep^2*m2 + ep*m1 + m0)*num(ep)*num(ep^3*n3 + ep^2*n2 + ep*n1 + n0);
Local test5 = + den(ep^3*m3 + ep^2*m2 + ep*m1 + m0)*num(ep^6*n6 + ep^5*n5 + ep^4*n4 + ep^3*n3 + ep^2*n2 + ep*n1 + n0);

FactArg num;
FactArg den;

SplitArg ((ep)) num;
Identify num(ep) = ep;
Transform num addargs(2,last);
FactArg num;

Print +s;
.end

507-3.frm

Symbol n6,...,n1;
*Symbol n1,...,n6;
Symbol ep;
CFunction num;

Local test4 =
    + num(ep*n1 + ep^2*n2 + ep^3*n3 + ep^4*n4 + ep^5*n5 + ep^6*n6)
   ;

FactArg num;

Print +s;
.end

507-4.frm

Symbol n6,...,n1;
*Symbol n1,...,n6;
Symbol ep;
CFunction num;

Local test4 =
    + num(ep*n1 + ep^2*n2 + ep^3*n3 + ep^4*n4 + ep^5*n5 + ep^6*n6)
   ;

FactArg num;
Transform num mulargs(1,last);
.sort

Print +s;
.end

507-5.frm

CF f;
S a,b,x;
L F = f(a*x^3 + 1) * f(x^2 + a*x + 1) * f(a*x^3 + b*x^2 + a*x + 1);
factarg f;
.end