[BACK]Return to de.rr CVS log [TXT][DIR] Up to [local] / OpenXM / src / asir-contrib / testing / noro

Annotation of OpenXM/src/asir-contrib/testing/noro/de.rr, Revision 1.1

1.1     ! noro        1: module de;
        !             2:
        !             3: localf split_lexgb;
        !             4: localf sort_lex_dec,sort_lex_inc;
        !             5: localf inverse_or_split, linear_dim;
        !             6: localf dp_monic_mod,monic_gb;
        !             7: localf membership_test;
        !             8: localf dp_chrem,intdptoratdp,intdpltoratdpl;
        !             9: localf comp_by_ht,dp_gr_mod,gr_chrem;
        !            10: localf construct_sqfrbasis;
        !            11:
        !            12: /*
        !            13:  * G : a 0-dim lex gb, reduced
        !            14:  * if there exists a non-monic g = g(x')x^k+...
        !            15:  * then g(x') is a zero divisor and Id(G) splits
        !            16:  */
        !            17:
        !            18: def split_lexgb(G,V)
        !            19: {
        !            20:        G = map(ptozp,G);
        !            21:        G = sort_lex_dec(G,V);
        !            22:        TG = G;
        !            23:        for ( ; TG != []; TG = cdr(TG) ) {
        !            24:                F = car(TG);
        !            25:                for ( TV = V; !deg(F,car(TV)); TV = cdr(TV) );
        !            26:                VF = car(TV);
        !            27:                DF = deg(F,VF);
        !            28:                CF = coef(F,DF,VF);
        !            29:                if ( type(CF) == 2 ) {
        !            30:                        L = inverse_or_split(V,G,CF);
        !            31:                        R = map(split_lexgb,L,V);
        !            32:                        return append(R[0],R[1]);
        !            33:                }
        !            34:        }
        !            35:        return [G];
        !            36: }
        !            37:
        !            38: /*
        !            39:  * V : a variable list
        !            40:  * Id : a 0-dim radical ideal represented by a lex basis
        !            41:  * F : a poly
        !            42:  * if F is a unit of Q[V]/Id, then 1/F is returned
        !            43:  * else F is a zero divisor and Id = (Id:F)\cap(Id+<F>)
        !            44:  *      [GB(Id:F),GB(Id+<F>)] is returned.
        !            45:  */
        !            46:
        !            47: def inverse_or_split(V,Id,F)
        !            48: {
        !            49:        Id = map(ptozp,Id);
        !            50:        N = length(V);
        !            51:        dp_ord(2);
        !            52:        set_field(Id,V,2);
        !            53:        DF = dptodalg(dp_ptod(F,V));
        !            54:        Ret = inv_or_split_dalg(DF);
        !            55:        if ( type(Ret) == 1 ) {
        !            56:                /* Ret = 1/F */
        !            57:                Dp = dalgtodp(Ret);
        !            58:                return dp_dtop(Dp[0],V)/Dp[1];
        !            59:        } else {
        !            60:                /* Ret = GB(Id:F) */
        !            61:                /* compute GB(Id+<f>) */
        !            62:                Gquo = append(map(ptozp,map(dp_dtop,Ret,V)),Id);
        !            63:                /* inter-reduction */
        !            64:                Gquo = nd_gr_postproc(Gquo,V,0,2,0);
        !            65:                DTotal = linear_dim(Id,V,2);
        !            66:                Dquo = linear_dim(Gquo,V,2);
        !            67:                Drem = DTotal-Dquo;
        !            68:                B = cons(F,Id);
        !            69:                Grem = gr_chrem(B,V,2|dim=Drem);
        !            70:                return [map(ptozp,Gquo),map(ptozp,Grem)];
        !            71:        }
        !            72: }
        !            73:
        !            74: /* add F(X,V) to Id(B) */
        !            75: /* returns a list of splitted ideals */
        !            76: /* B should be a triangular basis */
        !            77:
        !            78: def construct_sqfrbasis(F,X,B,V)
        !            79: {
        !            80:        if ( type(F) == 1 )
        !            81:                return [];
        !            82:        B = sort_lex_dec(B,V);
        !            83:        V1 = cons(X,V);
        !            84:        F = nd_nf(F,reverse(B),cons(X,V),2,0);
        !            85:        D = deg(F,X);
        !            86:        H = coef(F,D,X);
        !            87:        if ( type(H) == 2 ) {
        !            88:                Ret = inverse_or_split(V,B,H);
        !            89:                if ( type(Ret) == 4 ) {
        !            90:                        /* H != 0 on Id_nz, H = 0 on Id_z */
        !            91:                        B0=construct_sqfrbasis(F,X,Ret[0],V);
        !            92:                        B1=construct_sqfrbasis(F,X,Ret[1],V);
        !            93:                        return append(B0,B1);
        !            94:                } else
        !            95:                        F = nd_nf(F*Ret,reverse(B),cons(X,V),2,0);
        !            96:        }
        !            97:        B1 = cons(F,B);
        !            98:        /* F is monic */
        !            99:        M = minipoly(B1,V1,2,X,zzz);
        !           100:        S = sqfr(M); S = cdr(S);
        !           101:        if ( length(S) == 1 && car(S)[1] == 1 )
        !           102:                return [cons(F,B)];
        !           103:        else {
        !           104:                R = [];
        !           105:                for ( T = S; T != []; T = cdr(T) ) {
        !           106:                        G = nd_gr_trace(cons(subst(car(T)[0],zzz,X),B1),V1,1,1,2);
        !           107:                        R1 = split_lexgb(G,V1);
        !           108:                        R = append(R1,R);
        !           109:                }
        !           110:                return R;
        !           111:        }
        !           112: }
        !           113:
        !           114: def sort_lex_dec(B,V)
        !           115: {
        !           116:        dp_ord(2);
        !           117:        B = map(dp_ptod,B,V);
        !           118:        B = vtol(qsort(ltov(B),comp_by_ht));
        !           119:        B = map(dp_dtop,B,V);
        !           120:        return reverse(B);
        !           121: }
        !           122:
        !           123: def sort_lex_inc(B,V)
        !           124: {
        !           125:        dp_ord(2);
        !           126:        B = map(dp_ptod,B,V);
        !           127:        B = vtol(qsort(ltov(B),comp_by_ht));
        !           128:        B = map(dp_dtop,B,V);
        !           129:        return B;
        !           130: }
        !           131:
        !           132: def linear_dim(G,V,Ord)
        !           133: {
        !           134:        dp_ord(Ord);
        !           135:        MB = dp_mbase(map(dp_ptod,G,V));
        !           136:        return length(MB);
        !           137: }
        !           138:
        !           139: def membership_test(B,G,V,O)
        !           140: {
        !           141:        B = map(ptozp,B);
        !           142:        G = map(ptozp,G);
        !           143:        for ( T = B; T != []; T = cdr(T) )
        !           144:                if ( nd_nf(car(T),G,V,O,0) ) return 0;
        !           145:        return 1;
        !           146: }
        !           147:
        !           148: def gr_chrem(B,V,O)
        !           149: {
        !           150:        /* linear dimension */
        !           151:        Dim = getopt(dim);
        !           152:        if ( type(Dim) == -1 ) Dim = -1;
        !           153:
        !           154:        B = map(ptozp,B);
        !           155:        G = []; HS = []; Mod = 1;
        !           156:        for ( I = 0; ; I++ ) {
        !           157:                P = lprime(I);
        !           158:                GM = nd_gr(B,V,P,O);
        !           159:                if ( Dim >= 0 )
        !           160:                        if ( linear_dim(GM,V,O) != Dim ) continue;
        !           161:                L = monic_gb(GM,V,O,P); GM = L[0]; HSM = L[1];
        !           162:                G = dp_chrem(G,HS,Mod,GM,HSM,P);
        !           163:                Mod *= P;
        !           164:                if ( G != [] )
        !           165:                        HS = HSM;
        !           166:                M = idiv(isqrt(2*Mod),2);
        !           167:                R1 = intdpltoratdpl(G,Mod,M);
        !           168:                if ( R1 ) {
        !           169:                        if ( Found && R == R1
        !           170:                                && (GB=nd_gr_postproc(map(dp_dtop,R,V),V,0,O,1))
        !           171:                                && membership_test(B,GB,V,O) )
        !           172:                                break;
        !           173:                        else {
        !           174:                                R = R1; Found = 1;
        !           175:                        }
        !           176:                }
        !           177:        }
        !           178:        return GB;
        !           179: }
        !           180:
        !           181: def comp_by_ht(A,B)
        !           182: {
        !           183:        HA = dp_ht(A); HB = dp_ht(B);
        !           184:        if ( HA > HB )
        !           185:                return 1;
        !           186:        else if ( HA < HB )
        !           187:                return -1;
        !           188:        else
        !           189:                return 0;
        !           190: }
        !           191:
        !           192: def intdpltoratdpl(G,Mod,M)
        !           193: {
        !           194:        for ( R = []; G != []; G = cdr(G) ) {
        !           195:                T = intdptoratdp(car(G),Mod,M);
        !           196:                if ( !T )
        !           197:                        return 0;
        !           198:                else
        !           199:                        R = cons(T,R);
        !           200:        }
        !           201:        R = reverse(R);
        !           202:        return vtol(qsort(newvect(length(R),R),comp_by_ht));
        !           203: }
        !           204:
        !           205: def intdptoratdp(F,Mod,M)
        !           206: {
        !           207:        for ( T = F, N = 0; T; T = dp_rest(T), N++ );
        !           208:        C = newvect(N);
        !           209:        for ( I = 0, T = F; I < N; T = dp_rest(T), I++ ) {
        !           210:                L = inttorat(dp_hc(T),Mod,M);
        !           211:                if ( !L )
        !           212:                        return 0;
        !           213:                else
        !           214:                        C[I] = (L[0]/L[1])*dp_ht(T);
        !           215:        }
        !           216:        for ( R = 0, I = N-1; I >= 0; I-- )
        !           217:                R += C[I];
        !           218:        return R;
        !           219: }
        !           220:
        !           221: def dp_chrem(G,HS,Mod,GM,HSM,P)
        !           222: {
        !           223:        if ( G == [] )
        !           224:                return GM;
        !           225:        if ( HS != HSM )
        !           226:                return [];
        !           227:        R = [];
        !           228:        M1 = inv(Mod,P);
        !           229:        ModM1 = Mod*M1;
        !           230:        for ( ; G != []; G = cdr(G), GM = cdr(GM) ) {
        !           231:                E = car(G); EM = car(GM);
        !           232:                E1 = E+(EM-E)*ModM1;
        !           233:                R = cons(E1,R);
        !           234:        }
        !           235:        return reverse(R);
        !           236: }
        !           237:
        !           238: def monic_gb(G,V,O,P)
        !           239: {
        !           240:        dp_ord(O); setmod(P);
        !           241:        D = map(dp_ptod,G,V);
        !           242:        D = map(dp_monic_mod,D,P);
        !           243:        D = vtol(qsort(newvect(length(D),D),comp_by_ht));
        !           244:        return [D,map(dp_ht,D)];
        !           245: }
        !           246:
        !           247: def dp_monic_mod(F,P)
        !           248: {
        !           249:        FP = dp_mod(F,P,[]);
        !           250:        return dp_rat(FP/dp_hc(FP));
        !           251: }
        !           252:
        !           253: endmodule;
        !           254: end$

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>