Return to solve CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / lib

Annotation of OpenXM_contrib2/asir2000/lib/solve, Revision 1.1

1.1     ! noro        1: /* $OpenXM: OpenXM/src/asir99/lib/solve,v 1.1.1.1 1999/11/10 08:12:30 noro Exp $ */
        !             2: def kenzan(El,Sl)
        !             3: {
        !             4:        for ( Tl = El; Tl != []; Tl = cdr(Tl) ) {
        !             5:                if ( substv(car(Tl),Sl) ) {
        !             6:                        print("kenzan : error");
        !             7:                        return 0;
        !             8:                }
        !             9:        }
        !            10:        print("kenzan : ok");
        !            11:        return 1;
        !            12: }
        !            13:
        !            14: def substv(P,Sl)
        !            15: {
        !            16:        for ( A = P; Sl != []; Sl = cdr(Sl) )
        !            17:                A = subst(A,car(car(Sl)),car(cdr(car(Sl))));
        !            18:        return A;
        !            19: }
        !            20:
        !            21: def co(X,V,D)
        !            22: {
        !            23:        for ( I = 0; I < D; I++ )
        !            24:                X = diff(X,V);
        !            25:        return sdiv(subst(X,V,0),fac(D));
        !            26: }
        !            27:
        !            28: def solve(El,Vl)
        !            29: /*
        !            30:  * El : list of linear forms
        !            31:  * Vl : list of variable
        !            32:  */
        !            33: {
        !            34:        N = length(El); M = length(Vl);
        !            35:        Mat = newmat(N,M+1);
        !            36:        W = newvect(M+1); Index = newvect(N); Vs = newvect(M);
        !            37:        for ( I = 0, Tl = Vl; I < M; Tl = cdr(Tl), I++ )
        !            38:                Vs[I] = car(Tl);
        !            39:        for ( I = 0, Tl = El; I < N; Tl = cdr(Tl), I++ ) {
        !            40:                ltov(car(Tl),Vl,W);
        !            41:                for ( J = 0; J <= M; J++ )
        !            42:                        Mat[I][J] = W[J];
        !            43:        }
        !            44:        Tl = solvemain(Mat,Index,N,M); L = car(Tl); D = car(cdr(Tl));
        !            45:        if ( L < 0 )
        !            46:                return [];
        !            47:        for ( I = L - 1, S = []; I >= 0; I-- ) {
        !            48:                for ( J = Index[I]+1, A = 0; J < M; J++ ) {
        !            49:                        A += Mat[I][J]*Vs[J];
        !            50:                }
        !            51:                S = cons([Vs[Index[I]],-red((A+Mat[I][M])/D)],S);
        !            52:        }
        !            53:        if ( kenzan(El,S) )
        !            54:                return S;
        !            55:        else
        !            56:                return [];
        !            57:        return S;
        !            58: }
        !            59:
        !            60: def solvemain(Mat,Index,N,M)
        !            61: /*
        !            62:  *     Mat : matrix of size Nx(M+1)
        !            63:  *     Index : vector of length N
        !            64:  */
        !            65: {
        !            66:        for ( J = 0, L = 0, D = 1; J < M; J++ ) {
        !            67:                for ( I = L; I < N && !Mat[I][J]; I++ );
        !            68:                if ( I == N )
        !            69:                        continue;
        !            70:                Index[L] = J;
        !            71:                for ( K = 0; K <= M; K++ ) {
        !            72:                        T = Mat[I][K]; Mat[I][K] = Mat[L][K]; Mat[L][K] = T;
        !            73:                }
        !            74:                for ( I = L + 1, V = Mat[L][J]; I < N; I++ )
        !            75:                        for ( K = J, U = Mat[I][J]; K <= M; K++ )
        !            76:                                Mat[I][K] = sdiv(Mat[I][K]*V-Mat[L][K]*U,D);
        !            77:                D = V; L++;
        !            78:        }
        !            79:        for ( I = L; I < N; I++ )
        !            80:                for ( J = 0; J <= M; J++ )
        !            81:                        if ( Mat[I][J] )
        !            82:                                return -1;
        !            83:        for ( I = L - 2, W = newvect(M+1); I >= 0; I-- ) {
        !            84:                for ( J = 0; J <= M; J++ )
        !            85:                        W[J] = 0;
        !            86:                for ( G = I + 1; G < L; G++ )
        !            87:                        for ( H = Index[G], U = Mat[I][H]; H <= M; H++ )
        !            88:                                W[H] += Mat[G][H]*U;
        !            89:                for ( J = Index[I], U = Mat[I][J]; J <= M; J++ )
        !            90:                        Mat[I][J] = sdiv(Mat[I][J]*D-W[J],U);
        !            91:        }
        !            92:        return [L,D];
        !            93: }
        !            94:
        !            95: def length(L)
        !            96: {
        !            97:        for ( I = 0; L != []; L = cdr(L), I++ );
        !            98:        return I;
        !            99: }
        !           100:
        !           101: def ltov(P,VL,W)
        !           102: {
        !           103:        for ( I = 0, L = VL; L != []; L = cdr(L), I++ ) {
        !           104:                W[I] = co(P,car(L),1); P -= W[I]*car(L);
        !           105:        }
        !           106:        W[I] = P;
        !           107: }
        !           108: end$