#-
*	Program to compute the complete 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.
*	Don't forget to disable the Modulus in form.set!!!!!!
*
*	We need the mzv tables in separate files.
*	We assume that these tables only contain finite elements and that
*	duality has been applied. This is a bit slower but saves much space.
*
*	Program is an adaptation of hex.frm.
*	Made by J.Vermaseren, 4-aug-2008
*
*--#[ Settings :
*
#ifndef `WEIGHT'
#define WEIGHT "16"
#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 7
#case 8
#define GROUPING "8"
#break
#case 9
#case 10
#define GROUPING "16"
#break
#case 11
#case 12
#define GROUPING "32"
#break
#case 13
#case 14
#define GROUPING "64"
#break
#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
#default
#define GROUPING "512"
#break
#endswitch
#endif
*
#include hexc.h
#message Run with Weight = `WEIGHT', GROUPING = `GROUPING', ORDER = `ORDER', Rational calculus
#call readmzvs({`WEIGHT'-1})
.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.
*
Off Parallel;
L   FF = E(`WEIGHT'-2,1);
repeat id E(n?pos_,?a) = E(n-1,0,?a)+E(n-1,1,?a);
.sort
On Parallel;
#call duality(E,0)
.sort
DropCoefficient;
id  E(?a)=E(?a)*H(?a);
#call frombasis(`WEIGHT')
*
*	The following enforces a better ordering inside the equations.
*
#call convert(E,`WEIGHT',0)
#call convert(H,`WEIGHT',0)
*
*   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,`WEIGHT'/2
#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;
#do relation = Gen
.sort
InParallel;
Drop 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);
*
*   Empirical restrictions
*
id  E(1,x?,?a) = 0;
if ( count(fun0,1) == 2 );
elseif ( match(fun0(1)) );
elseif ( match(fun1(0)) );
elseif ( match(fun1(1)*fun0(2)) );
else;
    Discard;
endif;
*
*	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;
#do relation = Gen
.sort
InParallel;
Drop 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)
id  H(?a) = HH(?a);
id  E(?a) = H(?a);
B   H;
Print +f;
*--#] Output :
.sort
Off Parallel;
Global mzv`WEIGHT'c = FF;
if ( expression(FF) );
	ArgImplode,H;
	id	H(x1?,...,x{`WEIGHT'/3}?,x0?,?a) = 0;
	if ( count(HH,1) == 0 ) Discard;
endif;
B	H;
Print +f FF;
.store
save mzv`WEIGHT'.sav;
.end