Return to solve CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / lib

Annotation of OpenXM_contrib2/asir2000/lib/solve, Revision 1.1.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$