Return to sym CVS log

Up to [local] / OpenXM_contrib2 / asir2018 / lib

Annotation of OpenXM_contrib2/asir2018/lib/sym, Revision 1.1

1.1     ! noro        1: /* $OpenXM$ */
        !             2:
        !             3: /* functions for operations on symmetric polynomials */
        !             4:
        !             5: def gen_esym(N,V)
        !             6: {
        !             7:        R = newvect(N);
        !             8:        for ( I = 0, T = <<0>>; I < N; I++ )
        !             9:                T *= <<1>>+V[I]*<<0>>;
        !            10:        T = dp_rest(T);
        !            11:        for ( I = 0; I < N; I++, T = dp_rest(T) )
        !            12:                R[I] = dp_hc(T);
        !            13:        return R;
        !            14: }
        !            15:
        !            16: /*
        !            17:  This function expresses a symmetric polynomial F in terms of
        !            18:  the elementary symmetric polynomials.
        !            19:  If F is not symmetric it calls error().
        !            20:  It uses variables s1,s2,...,sn to represent elementary functions by default.
        !            21:  option :
        !            22:    vars=[x1,x2,...] => subset of vars(F).
        !            23:      F is assumed with respect to the specified variable set.
        !            24:    svars=[t1,t2,...] => #svars must be equal to #vars.
        !            25:      The specified variables are used instead of s1,s2,...,sn.
        !            26:  Example:
        !            27: [...] symmetric_reduction(x*y^2+x^2*y+y*z^2+y^2*z+z*x^2+z^2*x);
        !            28: -3*s3+s1*s2
        !            29: [...] symmetric_reduction(a*(x^3+y^3)+b*(x^2+y^2+x*y)|vars=[x,y]);
        !            30: (-3*s1*s2+s1^3)*a+(-s2+s1^2)*b
        !            31:  symmetric_reduction(s1^3+s2^3+s3^3|svars=[t1,t2,t3]);
        !            32:  t1^3-3*t2*t1+3*t3
        !            33: */
        !            34:
        !            35: def symmetric_reduction(F)
        !            36: {
        !            37:        V = getopt(vars);
        !            38:        if ( type(V) == -1 )
        !            39:                V = vars(F);
        !            40:        N = length(V);
        !            41:        S = getopt(svars);
        !            42:        if ( type(S) == -1 ) {
        !            43:                S = [];
        !            44:                for ( I = N; I >= 1; I-- )
        !            45:                        S = cons(strtov("s"+rtostr(I)),S);
        !            46:        }
        !            47:        dp_ord(2);
        !            48:        DF = dp_ptod(F,V);
        !            49:        ES = gen_esym(N,V);
        !            50:        while ( DF ) {
        !            51:                HE = dp_etov(DF);
        !            52:                for ( I = 1; I < N; I++ )
        !            53:                        if ( HE[I] > HE[I-1] )
        !            54:                                break;
        !            55:                if ( I < N )
        !            56:                        error("symmetric_reduction : input is not a symmetric polynomial");
        !            57:                K = N-1;
        !            58:                T = 1;
        !            59:                TT = 1;
        !            60:                while ( 1 ) {
        !            61:                        for ( ; K >= 0 && HE[K] == 0; K-- );
        !            62:                        if ( K < 0 )
        !            63:                                break;
        !            64:                        T *= ES[K]^HE[K];
        !            65:                        TT *= S[K]^HE[K];
        !            66:                        for ( I = 0; I < K; I++ )
        !            67:                                HE[I] -= HE[K];
        !            68:                        HE[K] = 0;
        !            69:                }
        !            70:                R += dp_hc(DF)*TT;
        !            71:                DF += -dp_hc(DF)*dp_ptod(T,V);
        !            72:        }
        !            73:        return R;
        !            74: }
        !            75:
        !            76: /* sum of all monomials of degree K with N variables */
        !            77: def monomials_of_eq_deg(N,K)
        !            78: {
        !            79:        W = newvect(N);
        !            80:        W[N-1] = K;
        !            81:        R = dp_vtoe(W);
        !            82:        if ( !K )
        !            83:                return R;
        !            84:        while ( next_exponent(W,K,N) )
        !            85:                R += dp_vtoe(W);
        !            86:        return R;
        !            87: }
        !            88:
        !            89: def next_exponent(W,K,N)
        !            90: {
        !            91:        for ( I = N-1; I >= 0 && W[I] == 0; I-- );
        !            92:        if ( I == 0 )
        !            93:                return 0;
        !            94:        else {
        !            95:                W[I-1]++;
        !            96:                T = W[I];
        !            97:                W[I] = 0;
        !            98:                W[N-1]=T-1;
        !            99:        }
        !           100:        return 1;
        !           101: }
        !           102:
        !           103: def sym_gb(V)
        !           104: {
        !           105:        N = length(V);
        !           106:        S = [];
        !           107:        for ( I = N; I >= 1; I-- )
        !           108:                S = cons(strtov("s"+rtostr(I)),S);
        !           109:        R = [];
        !           110:        for ( K = 1; K <= N; K++, V = cdr(V) ) {
        !           111:                T = monomials_of_eq_deg(N-K+1,K);
        !           112:                for ( J = 1; J <= K; J++ )
        !           113:                        T += (-1)^J*monomials_of_eq_deg(N-K+1,K-J)*S[J-1];
        !           114:                R = cons(dp_dtop(T,V),R);
        !           115:        }
        !           116:        return R;
        !           117: }
        !           118: end$