Return to sturm CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / lib

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

1.1     ! noro        1: /* $OpenXM: OpenXM/src/asir99/lib/sturm,v 1.1.1.1 1999/11/10 08:12:31 noro Exp $ */
        !             2: /* find intervals which include roots of a polynomial */
        !             3:
        !             4: #include "defs.h"
        !             5:
        !             6: def slice(P,XR,YR,WH) {
        !             7:        X = FIRST(XR); XMIN = SECOND(XR); XMAX = THIRD(XR);
        !             8:        Y = FIRST(YR); YMIN = SECOND(YR); YMAX = THIRD(YR);
        !             9:        W = FIRST(WH); H = SECOND(WH);
        !            10:        XS = (XMAX-XMIN)/W; YS = (YMAX-YMIN)/H;
        !            11:        T = ptozp(subst(P,X,X*XS+XMIN,Y,Y*YS+YMIN));
        !            12:        R = newvect(W+1);
        !            13:        for ( I = 0; I <= W; I++ ) {
        !            14:                S = sturm(subst(T,X,I));
        !            15:                R[I] = numch(S,Y,0)-numch(S,Y,H);
        !            16:        }
        !            17:        return R;
        !            18: }
        !            19:
        !            20: def slice1(P,XR,YR,WH) {
        !            21:        X = FIRST(XR); XMIN = SECOND(XR); XMAX = THIRD(XR);
        !            22:        Y = FIRST(YR); YMIN = SECOND(YR); YMAX = THIRD(YR);
        !            23:        W = FIRST(WH); H = SECOND(WH);
        !            24:        XS = (XMAX-XMIN)/W; YS = (YMAX-YMIN)/H;
        !            25:        T = ptozp(subst(P,X,X*XS+XMIN,Y,Y*YS+YMIN));
        !            26:        R = newvect(W+1);
        !            27:        for ( I = 0; I <= W; I++ ) {
        !            28:                S = sturm(subst(T,X,I));
        !            29:                R[I] = newvect(H+1);
        !            30:                seproot(S,Y,0,H,R[I]);
        !            31:        }
        !            32:        return R;
        !            33: }
        !            34:
        !            35: def seproot(S,X,MI,MA,R) {
        !            36:        N = car(size(S));
        !            37:        for ( I = MI; I <= MA; I++ )
        !            38:                if ( !(T = subst(S[0],X,I)) )
        !            39:                        R[I] = -1;
        !            40:                else
        !            41:                        break;
        !            42:        if ( I > MA )
        !            43:                return;
        !            44:        for ( J = MA; J >= MI; J-- )
        !            45:                if ( !(T = subst(S[0],X,J)) )
        !            46:                        R[J] = -1;
        !            47:                else
        !            48:                        break;
        !            49:        R[I] = numch(S,X,I); R[J] = numch(S,X,J);
        !            50:        if ( J <= I+1  )
        !            51:                return;
        !            52:        if ( R[I] == R[J] ) {
        !            53:                for ( K = I + 1; K < J; K++ )
        !            54:                        R[K] = R[I];
        !            55:                return;
        !            56:        }
        !            57:        T = idiv(I+J,2);
        !            58:        seproot(S,X,I,T,R);
        !            59:        seproot(S,X,T,J,R);
        !            60: }
        !            61:
        !            62: /* compute the sturm sequence of P */
        !            63:
        !            64: def sturm(P) {
        !            65:        V = var(P); N = deg(P,V); T = newvect(N+1);
        !            66:        G1 = T[0] = P; G2 = T[1] = ptozp(diff(P,var(P)));
        !            67:        for ( I = 1, H = 1, X = 1; ; ) {
        !            68:                if ( !deg(G2,V) )
        !            69:                        break;
        !            70:                D = deg(G1,V)-deg(G2,V);
        !            71:                if ( (L = LCOEF(G2)) < 0 )
        !            72:                        L = -L;
        !            73:                if ( !(R = -srem(L^(D+1)*G1,G2)) )
        !            74:                        break;
        !            75:                if ( type(R) == 1 ) {
        !            76:                        T[++I] = (R>0?1:-1); break;
        !            77:                }
        !            78:                M = H^D; G1 = G2;
        !            79:                G2 = T[++I] = sdiv(R,M*X);
        !            80:                if ( (X = LCOEF(G1)) < 0 )
        !            81:                        X = -X;
        !            82:                H = X^D*H/M;
        !            83:        }
        !            84:        S = newvect(I+1);
        !            85:        for ( J = 0; J <= I; J++ )
        !            86:                S[J] = T[J];
        !            87:        return S;
        !            88: }
        !            89:
        !            90: def numch(S,V,A) {
        !            91:        N = car(size(S));
        !            92:        for ( T = subst(S[0],V,A), I = 1, C = 0; I < N; I++ ) {
        !            93:                T1 = subst(S[I],V,A);
        !            94:                if ( !T1 )
        !            95:                        continue;
        !            96:                if ( (T1 > 0 && T < 0) || (T1 < 0 && T > 0) )
        !            97:                        C++;
        !            98:                T = T1;
        !            99:        }
        !           100:        return C;
        !           101: }
        !           102:
        !           103: def numch0(S,V,A,T) {
        !           104:        N = car(size(S));
        !           105:        for ( I = 1, C = 0; I < N; I++ ) {
        !           106:                T1 = subst(S[I],V,A);
        !           107:                if ( !T1 )
        !           108:                        continue;
        !           109:                if ( (T1 > 0 && T < 0) || (T1 < 0 && T > 0) )
        !           110:                        C++;
        !           111:                T = T1;
        !           112:        }
        !           113:        return C;
        !           114: }
        !           115: end;