File: [local] / OpenXM_contrib2 / asir2000 / lib / sym (download)
Revision 1.1, Fri Mar 9 01:10:58 2001 UTC (23 years, 6 months ago) by noro
Branch: MAIN
CVS Tags: R_1_3_1-2, RELEASE_1_3_1_13b, RELEASE_1_2_3_12, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, KNOPPIX_2006, HEAD, DEB_REL_1_2_3-9
Added a library file 'sym', which contains several functions related to
symmetric polynomial operations.
|
/* $OpenXM: OpenXM_contrib2/asir2000/lib/sym,v 1.1 2001/03/09 01:10:58 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$