*
*   Master include file for the hec.frm program that computes H(R(...),1)
*   for a given weight WEIGHT.
*   File by J. Vermaseren 15-jul-2008
*
*--#[ clean :
*
#procedure clean(type)
*
#ifdef `WHERE'
#message doing `WHERE'
#endif
*
*   The handling of the equations once they have been properly constructed
*   First convert back to 'table notation' of equation 8 in the paper
*
  id  H(R(?a),1) = H(?a);
  id  HH(R(?a),1) = HH(?a);
  if ( count(H,1,ln2,1,Sinf,1) > 1 ) Multiply replace_(H,HH);
  id  HH(?a) = 0;
  Multiply replace_(HH,H);
  #call makeHfinite
  .sort: after makeHfinite(`DEPTH'-`type');
  InParallel;
#if ( {2*`DEPTH'} == `WEIGHT' )
  #do k = 1,`WEIGHT'-1
    id  H(n1?,...,n`k'?) = 0;
  #enddo
  #call duality(H,0)
  .sort: finite(`DEPTH'-`type');
  InParallel;
#endif
*
*   Now some trick to be able to substitute the H-functions of the proper
*   weight by the contents of the appropriate bracket in FF.
*
*  repeat id  H(?a,0,?b) = H(?a,3,?b);
#call convert(H,`WEIGHT',0)
  id  H(?a) = R(E(?a));
  id  R(n?$n) = FF[$n];
*
*   Bracket in H so that we may take the first bracket
*
  B   H;
ModuleOption,local,$n;
 .sort: setup(`DEPTH'-`type')-`$count';
ClearTable Hfill;
Off Statistics;
Off Parallel;
*
*	Now we have `$jj' equations
*
  #$v = 0;
  #do k1 = 1,`$jj'
  #if ( termsin(F`k1') != 0 )
	InParallel F1,...,F`$jj';
    Skip F`k1';
	#$v = $v + 1;
    #$vna`$v' = `k1';
    #$fb`$v' = FirstBracket_(F`k1');
    #$fv = F`k1'[$fb`$v'];
	#$fv = 1/(`$fv');
    #$fr`$v' = -F`k1'*(`$fv')+`$fb`$v'';
    id `$fb`$v'' = `$fr`$v'';
	B	H;
	.sort:Preparing `$dcount';
    #$dcount = $dcount-1;
  #endif
  #enddo
  #if ( `$v' > 0 )
	InParallel;
	id	H(?a) = Hfill(?a);
	B	Hfill;
	.sort:Setting Hfill `$dcount';
    #do k1 = 1,`$v'
      #$fb`k1' = FirstBracket_(F`$vna`k1'');
      #$fv = F`$vna`k1''[$fb`k1'];
	  #$fv = 1/(`$fv');
      #$fr`k1' = (-F`$vna`k1''*(`$fv')+`$fb`k1'')*replace_(Hfill,H);
    #enddo
    #do vv = 1,`$v'
      Fill,`$fb`vv'' = `$fr`vv'';
    #enddo
On Parallel;
    On Statistics;
    Drop;
    Ndrop FF;
    Unhide FF;
    id	H(?a) = Hfill(?a);
	id	Hfill(?a) = H(?a);
    B+   E;
    .sort:substitution(`DEPTH'-`type')-{`$dcount'+1};
Off Statistics;
*Off Parallel;
	InParallel;
	Hide FF;
  #else
Off Statistics;
*Off Parallel;
	InParallel;
    Drop;
    Ndrop FF;
  #endif
#endprocedure
*
*--#] clean :
*--#[ solveshuffle :
*
#procedure solveshuffle
*
*   Procedure for equations based on the shuffle algebra.
*   Based on the reduction to the basis (see appendix A).
*   EE*EE defines an equation. (HH*HH-H*H = 0)
*   clean does the rest (see in clean (above))
*
*Print +f "<sh> %t";
id,many,fun1?funs(?a) = 1;
  .sort
  InParallel;
  id  E(?a)*E(?b) = HH(?a)*HH(?b) - H(?a)*H(?b);
  id  HH(?a) = 0;
  Shuffle,H;
  .sort:Shuffling-`$count';
  InParallel;
  #call clean(sh)
#endprocedure
*
*--#] solveshuffle : 
*--#[ solvestuffle :
*
#procedure solvestuffle
*
*	Procedure for equations for H in 1. We convert to Z sums and
*	apply the basisz routine there.
*   Then clean does the rest (see in clean (above))
*
id,many,fun1?funs(?a) = 1;
.sort
  InParallel;
  id  E(?a)*E(?b) = HH(?a)*HH(?b) - H(?a)*H(?b);
  id  HH(?a) = 0;
#ifdef `GENERICSTUFFLE'
*
  ArgImplode H;
*
*	Convert to Z
*
  repeat;
	id,once,H(?a) = R1*R2(?a);
	repeat id R1(?a)*R2(?b,n1?,n2?) = sig_(n2)*R1(sig_(n1)*n2,?a)*R2(?b,n1);
	id	R1(?a)*R2(n1?) = sig_(n1)*Z(n1,?a);
  endrepeat;
*
  Stuffle,Z+;
*
.sort:Stuffling-`$count';
  InParallel;
*
*	Convert back to H
*
  repeat;
	id,once,Z(?a) = R1*R2(?a)*R3(1);
	repeat id R1(?a)*R2(n1?,?b)*R3(x?) = sig_(n1)*x*R1(?a,n1*x)*R2(?b)*R3(x*sig_(n1));
	id	R1(?a)*R2*R3(x?) = H(?a);
  endrepeat;
*
  ArgExplode H;
#else
*
*	Stuffle for when all indices are positive.
*
	ArgImplode H;
	Stuffle,H+;
	.sort:Stuffling-`$count';
	InParallel;
	ArgExplode H;
#endif
*
  #call clean(st)
#endprocedure
*
*--#] solvestuffle : 
*--#[ duality :
*
#procedure duality(H,par)
*
*	Applies duality on finite non-alternating H values in one
*	If par == 0 we take a WEIGHT specific (faster) algorithm
*	If par == 1 we take a general (slower) algorithm
*
#switch `par'
#case 0
  repeat;
   if ( count(`H',1) );
	id,once,`H'(n1?,...,n`WEIGHT'?) =
				<x1/x1^n`WEIGHT'>*...*<x`WEIGHT'/x`WEIGHT'^n1>
				*R3(n1,...,n`WEIGHT');
	if ( count(<x1,2>,...,<x`WEIGHT',2>) > `WEIGHT' );
		id,many,x?{x1,...,x`WEIGHT'} = 1;
		id	R3(?a) = R4(?a);
	elseif ( count(<x1,2>,...,<x`WEIGHT',2>) < `WEIGHT' );
		id,<x1^n1?>*...*<x`WEIGHT'^n`WEIGHT'?> = R4(n1,...,n`WEIGHT');
		id	R3(?a) = 1;
	else;
		id,<x1^n1?>*...*<x`WEIGHT'^n`WEIGHT'?> = R3(n1,...,n`WEIGHT');
		id,R3(?a)*R3(?b) = R4(?a);
	endif;
   endif;
  endrepeat;
  id  R4(?a) = `H'(?a);
#break
#case 1
	repeat;
	  if ( match(`H'(x1?,x2?,x3?,?a)) );
		id,once,`H'(x1?,x2?,x3?,?a) = R1(x1,x2,x3,?a)*R2*R3(x1,x2,x3,?a);
		repeat;
			id	R1(0,?a)*R2(?b) = R1(?a)*R2(1,?b)*xc;
			id	R1(1,?a)*R2(?b) = R1(?a)*R2(0,?b)/xc;
		endrepeat;
		id	R1 = 1;
		if ( count(xc,1) == 0 );
			id	R2(?a) = R3(?a);
			id	R3(?a)*R3(?b) = R4(?a);
		elseif ( count(xc,1) > 0 );
			id	R2(?a)*R3(?b) = R4(?b);
			id	xc^n? = 1;
		else;
			id	R2(?a)*R3(?b) = R4(?a);
			id	xc^n? = 1;
		endif;
	  endif;
	endrepeat;
	id	R4(?a) = `H'(?a);
#break
#endswitch
#endprocedure
*
*--#] duality : 
*--#[ makeHfinite :
*
#procedure makeHfinite
*
*	Extracts the divergent parts from a sum
*	We use single stuffles
*	Z(st(1,R(n,1)),m)+Z(R(n,1),m1,st(1,restofm)) = Z(1)*Z(R(n,1),m)
*	Z(R(n+1,1),m)*(n+1) + Z(stp(1,R(n,1),m)+Z(R(n,1),m1,st(1,restofm))
*			 = Z(1)*Z(R(n,1),m)
*
*	This routine constitues a major 'inefficiency'
*	We eliminate all powers of Sinf
*
*	id	H(1) = Sinf;
	id	H(1) = 0;
if ( match(H(1,?a)) );
	ArgImplode H;
	id	H(?a) = Z(?a);
endif;
#$jc = 0;
#do i = 1,1
#$jc = $jc+1;
	if ( match(Z(1,1,?a)) );
		id,once,Z(1,?a) = Za(1,?a);
		Multiply ZZZ;
		repeat id ZZZ(?a)*Za(1,?b) = ZZZ(?a,1)*Za(?b);
*		id	ZZZ(1,?a)*Za(?b) = (Sinf*Za(?a,?b)-Za(1)*Za(?a,?b)*ZZZ(1,?a))*modinverse(nargs_(?a)+1);
		id	ZZZ(1,?a)*Za(?b) = (-Za(1)*Za(?a,?b)*ZZZ(1,?a))*modinverse(nargs_(?a)+1);
		Stuffle,Za+;
		id	ZZZ(?a)*Za(?a,?b) = 0;
		id	ZZZ(?a) = 1;
		id	Za = 1;
		id	Za(?a) = Z(?a);
		if ( match(Z(1,?a)) ) redefine i "0";
	elseif ( match(Z(1,?a)) );
		id,once,Z(1,?a) = Za(1,?a)
				-Za(1)*Za(?a)
*				+Sinf*Za(?a)
					;
		Stuffle,Za+;
		if ( match(Z(1,?a)) ) redefine i "0";
		id	Za(?a) = Z(?a);
	endif;
.sort:makeHfinite(`$jc');
InParallel;
#enddo
*repeat;
id	Z(?a) = H(?a);
ArgExplode H;
#endprocedure
*
*--#] makeHfinite : 
*--#[ invtable :
*
#procedure invtable(modinverse,N)
*
*	Makes a table of the first N fractions 1/i.
*	This is for modular calculus. Hence we assume that the modulus statement
*	is active. It avoids having to use the extended Euclidean algorithm
*	over and over again.
*	We assume that there is a table declaration
*	Table modinverses(1:`N');
*
.sort
Skip;
L	Finvtable = sum_(i,1,`N',x^i);
id	x^j? = x^j/j;
B	x;
.sort
FillExpression `modinverse' = Finvtable(x);
Skip;
Drop Finvtable;
.sort
*
#endprocedure
*
*--#] invtable : 
*--#[ convert :
*
#procedure convert(H,w,par)
#if ( `par' == 0 )
	id `H'(x?,?a) = `H'(1-x,?a);
	repeat id `H'(x1?,x?,?a) = `H'(2*x1+1-x,?a);
#elseif ( `par' == 1 )
	#do i = 1,`w'-1
		id	`H'(x?,?a) = `H'((x-mod_(x,2))/2,1-mod_(x,2),?a);
	#enddo
	id `H'(x?,?a) = `H'(1-x,?a);
#endif
#endprocedure
*
*--#] convert :
*--#[ Startup :
*
CF	fun0,fun1,...,fun20,funa0,funa1,funa2; * Functions to help with ordering the equarion
Set funs:fun0,fun1,...,fun20,funa0,funa1,funa2;
#ifdef `MODULUS'
Modulus,plusmin,nodollars,`MODULUS';
#endif
AutoDeclare Symbols i,j,l,m,n,x,y,h,e,s,z;
AutoDeclare CF S,H,Z,K,sum,den,sign;
CF	R,RR,R1,R2,R3,R4,fff(s),ffs;
CF  E,EE,Eb,EN,binom,acc1,Markett,A,B;
CF	keeplyndonfun0,keeplyndonfun1,keeplyndonfun2;
S	n,k,a,b,c,aa,ab,ac,sa,sb,sc;
S	Sinf,ln2,z2;
Table,sparse,Hfill(1);
*Table,sparse,Hfill(`WEIGHT');
Table,modinverse(1:{`WEIGHT'*2});
*
#call invtable(modinverse,{`WEIGHT'*2})
*
On Parallel;
On HighFirst;
On FewerStatistics;
Format nospaces;
*
*--#] Startup :