#-
*	Program to compute a Lyndon basis for the nonalternating MZV's
*	of a given weight (indicated by the parameter WEIGHT).
*
*	Settings via form.set because they are rather detailed.
*	And we carry this program between computers of different types.
*
*	Program by J.Vermaseren, 15-jul-2008
*
*--#[ Settings :
*
#ifndef `WEIGHT'
#define WEIGHT "20"
#endif
#define MAXDEPTH "20"
#define MODULUS "2147479273"
#if ( {2*`MAXDEPTH'} > `WEIGHT' )
#redefine MAXDEPTH "{`WEIGHT'/2}"
#endif
#ifndef `ORDER'
#define ORDER "X1"
#endif
*
*   Below is the optimum as determined in a number of trial runs.
*   It isn't 100% clear though whether this optimum has survived
*   later changes in orderings.
*
#ifndef `GROUPING'
#switch `WEIGHT'
#case 13
#case 14
#define GROUPING "64"
#case 15
#case 16
#define GROUPING "128"
#break
#case 17
#case 18
#define GROUPING "256"
#break
#case 19
#case 20
#define GROUPING "512"
#break
#case 21
#case 22
#define GROUPING "1024"
#break
#case 23
#case 24
#define GROUPING "2048"
#break
#case 25
#case 26
#define GROUPING "4096"
#break
#case 27
#case 28
#case 29
#case 30
#define GROUPING "6144"
#break
#default
#define GROUPING "512"
#break
#endswitch
#endif
*
#include hex.h
#ifdef `MODULUS'
#message Run with Weight = `WEIGHT', GROUPING = `GROUPING', ORDER = `ORDER', Modulus = `MODULUS'
#else
#message Run with Weight = `WEIGHT', GROUPING = `GROUPING', ORDER = `ORDER', Rational calculus
#endif
.global
*--#] Settings :
*--#[ Master Equation :
*
*   Construct the expression FF that will contain the final result.
*	We only do the finite elements.
*   At first it will contain 2^{`WEIGHT'-2} objects but duality will
*	remove either half or almost half the elements.
*
S	xcount(:`MAXDEPTH');
Off Parallel;
L   FF = E(`WEIGHT'-2,1)*xcount;
repeat id E(n?pos_,?a) = E(n-1,0,?a)+E(n-1,1,?a)*xcount;
.sort
On Parallel;
#call duality(E,0)
id	xcount^n? = 1;
.sort
DropCoefficient;
#call convert(E,`WEIGHT',0)
id  E(?a)=E(?a)*H(?a);
*
*	The following enforces a better ordering inside the equations.
*
*repeat id   E(?a,0,?b) = E(?a,3,?b);
*repeat id   H(?a,0,?b) = H(?a,3,?b);
*
*   The bracket is essential!
*
B+   E;
.sort
Off Statistics;
Hide FF;
.sort
*
*   $count will tell what is the number of the identity we are constructing
*   $dcount tell how many H-function have not been eliminated yet.
*	$jj counts the elements in the group that we treat simultaneously.
*
#$count = 0;
#$dcount = termsin_(FF) - 1;
#$jj = 0;
*
*--#] Master Equation :
*--#[ The relations :
*
#do DEPTH = 2,`MAXDEPTH'
#message entering DEPTH = `DEPTH'
*
*--#[ Stuffles :
*
*#message doing stuffles (`STUFFLE')({`DEPTH'-`STUFFLE'})
#message doing stuffles
*
Off Parallel;
L	Gen = E(`DEPTH')*EE({`WEIGHT'-`DEPTH'-1});
repeat id E(x?{>1},?a) = E(x-1,1,?a);
repeat id EE(x?pos_,?a) = EE(x-1,0,?a);
id	EE(?a) = E(?a);
shuffle,E;
id	E(?a,0) = 0;
id	E(1,?a) = 0;
.sort
On Parallel;
ArgImplode,E;

id	E(?a) = E(?a)*fff(?a);
Multiply replace_(fff,ffs);
#do i = 1,`DEPTH'
id	ffs(x?,?a) = ffs(?a)*fun`i'(-x);
#enddo
id	ffs(?a) = 1;
*----------------------------------------------------------------
*id	E(x1?,...,x`DEPTH'?) =
*#do i = 1,`DEPTH'/2
*		+E(x1,...,x`i')*E(<x{`i'+1}>,...,<x`DEPTH'>)*funa0(`i')
*#enddo
*		;
*----------------------------------------------------------------
Multiply,(
#do i = 1,`DEPTH'/2
	+EE(`i')*funa0(`i')
#enddo
		);
repeat id EE(x?pos_,?a)*E(x1?,?b) = EE(x-1,?a,x1)*E(?b);
id	EE(0,?a) = E(?a);
*----------------------------------------------------------------
id	E(1,?a) = 0;
.sort
DropCoefficient;
id	E(x1?,?a)*E(x2?,?b) = E(x1,?a)*E(x2,?b)*funa2(-min_(x1,x2));
ArgExplode,E;
.sort
Off Parallel;
InParallel;
Drop Gen;
#do relation = Gen
	#$jj = $jj+1;
	L	F`$jj' = `relation';
	#if ( `$jj' >= `GROUPING' )
		#call solvestuffle
		#$jj = 0;
	#endif
#enddo
#call solvestuffle
#$jj = 0;
.sort
*
*--#] Stuffles :
*--#[ Shuffles :
*
#do SHUFFLE = 1,1
#message doing shuffles (`SHUFFLE')({`DEPTH'-`SHUFFLE'})
*
Off Parallel;
L	Gen = E(`DEPTH')*EE({`WEIGHT'-`DEPTH'-1});
repeat id E(x?{>1},?a) = E(x-1,1,?a);
repeat id EE(x?pos_,?a) = EE(x-1,0,?a);
id	EE(?a) = E(?a);
shuffle,E;
id	E(?a,0) = 0;
id	E(1,1,?a) = 0;
ArgImplode,E;
Multiply EE;
#do i = 1,`SHUFFLE'
	id	EE(?a)*E(x?,?b) = EE(?a,x)*E(?b);
#enddo
id	EE(?a) = E(?a);
id	E(1,1,?a) = 0;
id	EE(1,?a)*E(1,?b) = 0;
.sort
DropCoefficient;
*
*	Look for the number of ones. The fewer the better.
*
#switch `ORDER'
#case X1
#do i = 1,1
id	E(?a) = E(R,R(?a));
repeat;
	id E(R(?a),R(`i',?b)) = E(R(?a,`i'),R(?b))*x;
	id E(R(?a),R(x1?!{`i',`i'},?b)) = E(R(?a,x1),R(?b));
endrepeat;
id	E(R(?a),R) = E(?a);
id	x^n? = fun`i'(n);
#enddo
#break
#case X2
#do i = 1,`WEIGHT'/2
id	E(?a) = E(R,R(?a));
repeat;
	id E(R(?a),R(`i',?b)) = E(R(?a,`i'),R(?b))*x;
	id E(R(?a),R(x1?!{`i',`i'},?b)) = E(R(?a,x1),R(?b));
endrepeat;
id	E(R(?a),R) = E(?a);
id	x^n? = fun`i'(n);
#enddo
#break
#endswitch
*
*	The next sets the lower weight in the single depth object first.
*
id	E(x1?) = E(x1)*fun0(x1);
*
*	Look for sum of squares of indices
*
id	E(x1?,?a) = E(x1^2,R(x1),R(?a));
repeat id E(x?,R(?a),R(x1?,?b)) = E(x+x1^2,R(?a,x1),R(?b));
id	E(x?,R(?a),R) = E(x,?a);
id	E(x1?,?a)*E(x2?,?b) = E(?a)*E(?b)*funa0(-max_(x1,x2));
*
*	Look for the leading index.
*
id	E(x1?,?a)*E(x2?,?b) = E(x1,?a)*E(x2,?b)*funa1(-min_(x1,x2))*funa2(-x1-x2);
ArgExplode,E;
.sort
#$jj = 0;
InParallel;
Drop Gen;
#do relation = Gen
	#$jj = $jj+1;
	L	F`$jj' = `relation';
	#if ( `$jj' >= `GROUPING' )
		#call solveshuffle
		#$jj = 0;
	#endif
#enddo
#call solveshuffle
#$jj = 0;
.sort
#enddo
*
*--#] Shuffles :
*
#enddo
*
*--#] The relations :
*--#[ Output :
.sort
Drop;
ndrop FF;
.sort
On Statistics;
On Parallel;
Format nospaces;
Format 80;
UnHide FF;
*
*   In the output the bracket indicates the H-function and the contents of
*   the bracket what it is equal to.
*   The remaining H functions are renamed as HH. This way they can be
*   found easily in an editor.
*
#call convert(E,`WEIGHT',1)
#call convert(H,`WEIGHT',1)
*repeat id   E(?a,3,?b) = E(?a,0,?b);
*repeat id   H(?a,3,?b) = H(?a,0,?b);
id  H(?a) = HH(?a);
id  E(?a) = H(?a);
B   H;
Print +f;
*--#] Output :
.sort
Off Parallel;
B	HH;
Print[] +f;
.end