/* $OpenXM: OpenXM_contrib2/asir2018/lib/sym,v 1.1 2018/09/19 05:45:08 noro Exp $ */

/* functions for operations on symmetric polynomials */

def gen_esym(N,V)
{
	R = newvect(N);
	for ( I = 0, T = <<0>>; I < N; I++ )
		T *= <<1>>+V[I]*<<0>>;
	T = dp_rest(T);
	for ( I = 0; I < N; I++, T = dp_rest(T) )
		R[I] = dp_hc(T);
	return R;
}

/*
 This function expresses a symmetric polynomial F in terms of
 the elementary symmetric polynomials.
 If F is not symmetric it calls error().
 It uses variables s1,s2,...,sn to represent elementary functions by default.
 option :
   vars=[x1,x2,...] => subset of vars(F).
     F is assumed with respect to the specified variable set.
   svars=[t1,t2,...] => #svars must be equal to #vars.
     The specified variables are used instead of s1,s2,...,sn.
 Example:
[...] symmetric_reduction(x*y^2+x^2*y+y*z^2+y^2*z+z*x^2+z^2*x);
-3*s3+s1*s2
[...] symmetric_reduction(a*(x^3+y^3)+b*(x^2+y^2+x*y)|vars=[x,y]);
(-3*s1*s2+s1^3)*a+(-s2+s1^2)*b
 symmetric_reduction(s1^3+s2^3+s3^3|svars=[t1,t2,t3]);
 t1^3-3*t2*t1+3*t3
*/

def symmetric_reduction(F)
{
	V = getopt(vars);
	if ( type(V) == -1 )
		V = vars(F);
	N = length(V);
	S = getopt(svars);
	if ( type(S) == -1 ) {
		S = [];
		for ( I = N; I >= 1; I-- )
			S = cons(strtov("s"+rtostr(I)),S);
	}
	dp_ord(2);
	DF = dp_ptod(F,V);
	ES = gen_esym(N,V);	
	while ( DF ) {
		HE = dp_etov(DF);
		for ( I = 1; I < N; I++ )
			if ( HE[I] > HE[I-1] )
				break;
		if ( I < N )
			error("symmetric_reduction : input is not a symmetric polynomial");
		K = N-1;
		T = 1;
		TT = 1;
		while ( 1 ) {
			for ( ; K >= 0 && HE[K] == 0; K-- );
			if ( K < 0 )
				break;
			T *= ES[K]^HE[K];
			TT *= S[K]^HE[K];
			for ( I = 0; I < K; I++ )
				HE[I] -= HE[K];
			HE[K] = 0;
		}
		R += dp_hc(DF)*TT;
		DF += -dp_hc(DF)*dp_ptod(T,V);
	}
	return R;
}

/* sum of all monomials of degree K with N variables */
def monomials_of_eq_deg(N,K)
{
	W = newvect(N);
	W[N-1] = K;
	R = dp_vtoe(W);
	if ( !K )
		return R;
	while ( next_exponent(W,K,N) )
	  	R += dp_vtoe(W);
	return R;
}

def next_exponent(W,K,N)
{
	for ( I = N-1; I >= 0 && W[I] == 0; I-- );
	if ( I == 0 )
		return 0;
	else {
		W[I-1]++;
		T = W[I];
		W[I] = 0;
		W[N-1]=T-1;
	}
	return 1;
}

def sym_gb(V)
{
	N = length(V);
	S = [];
	for ( I = N; I >= 1; I-- )
		S = cons(strtov("s"+rtostr(I)),S);
	R = [];
	for ( K = 1; K <= N; K++, V = cdr(V) ) {
		T = monomials_of_eq_deg(N-K+1,K);
		for ( J = 1; J <= K; J++ )
			T += (-1)^J*monomials_of_eq_deg(N-K+1,K-J)*S[J-1];
		R = cons(dp_dtop(T,V),R);
	}
	return R;
}
end$