[BACK]Return to de CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2018 / lib

Annotation of OpenXM_contrib2/asir2018/lib/de, 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,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,Dim)
        !           149: {
        !           150:        B = map(ptozp,B);
        !           151:        G = []; HS = []; Mod = 1;
        !           152:        for ( I = 0; ; I++ ) {
        !           153:                P = lprime(I);
        !           154:                GM = nd_gr(B,V,P,O);
        !           155:                if ( linear_dim(GM,V,O) != Dim ) continue;
        !           156:                L = monic_gb(GM,V,O,P); GM = L[0]; HSM = L[1];
        !           157:                G = dp_chrem(G,HS,Mod,GM,HSM,P);
        !           158:                Mod *= P;
        !           159:                if ( G != [] )
        !           160:                        HS = HSM;
        !           161:                M = idiv(isqrt(2*Mod),2);
        !           162:                R1 = intdpltoratdpl(G,Mod,M);
        !           163:                if ( R1 ) {
        !           164:                        if ( Found && R == R1
        !           165:                                && (GB=nd_gr_postproc(map(dp_dtop,R,V),V,0,O,1))
        !           166:                                && membership_test(B,GB,V,O) )
        !           167:                                break;
        !           168:                        else {
        !           169:                                R = R1; Found = 1;
        !           170:                        }
        !           171:                }
        !           172:        }
        !           173:        return GB;
        !           174: }
        !           175:
        !           176: def comp_by_ht(A,B)
        !           177: {
        !           178:        HA = dp_ht(A); HB = dp_ht(B);
        !           179:        if ( HA > HB )
        !           180:                return 1;
        !           181:        else if ( HA < HB )
        !           182:                return -1;
        !           183:        else
        !           184:                return 0;
        !           185: }
        !           186:
        !           187: def intdpltoratdpl(G,Mod,M)
        !           188: {
        !           189:        for ( R = []; G != []; G = cdr(G) ) {
        !           190:                T = intdptoratdp(car(G),Mod,M);
        !           191:                if ( !T )
        !           192:                        return 0;
        !           193:                else
        !           194:                        R = cons(T,R);
        !           195:        }
        !           196:        R = reverse(R);
        !           197:        return vtol(qsort(newvect(length(R),R),comp_by_ht));
        !           198: }
        !           199:
        !           200: def intdptoratdp(F,Mod,M)
        !           201: {
        !           202:        for ( T = F, N = 0; T; T = dp_rest(T), N++ );
        !           203:        C = newvect(N);
        !           204:        for ( I = 0, T = F; I < N; T = dp_rest(T), I++ ) {
        !           205:                L = inttorat(dp_hc(T),Mod,M);
        !           206:                if ( !L )
        !           207:                        return 0;
        !           208:                else
        !           209:                        C[I] = (L[0]/L[1])*dp_ht(T);
        !           210:        }
        !           211:        for ( R = 0, I = N-1; I >= 0; I-- )
        !           212:                R += C[I];
        !           213:        return R;
        !           214: }
        !           215:
        !           216: def dp_chrem(G,HS,Mod,GM,HSM,P)
        !           217: {
        !           218:        if ( G == [] )
        !           219:                return GM;
        !           220:        if ( HS != HSM )
        !           221:                return [];
        !           222:        R = [];
        !           223:        M1 = inv(Mod,P);
        !           224:        ModM1 = Mod*M1;
        !           225:        for ( ; G != []; G = cdr(G), GM = cdr(GM) ) {
        !           226:                E = car(G); EM = car(GM);
        !           227:                E1 = E+(EM-E)*ModM1;
        !           228:                R = cons(E1,R);
        !           229:        }
        !           230:        return reverse(R);
        !           231: }
        !           232:
        !           233: def monic_gb(G,V,O,P)
        !           234: {
        !           235:        dp_ord(O); setmod(P);
        !           236:        D = map(dp_ptod,G,V);
        !           237:        D = map(dp_monic_mod,D,P);
        !           238:        D = vtol(qsort(newvect(length(D),D),comp_by_ht));
        !           239:        return [D,map(dp_ht,D)];
        !           240: }
        !           241:
        !           242: def dp_monic_mod(F,P)
        !           243: {
        !           244:        FP = dp_mod(F,P,[]);
        !           245:        return dp_rat(FP/dp_hc(FP));
        !           246: }
        !           247:
        !           248: endmodule;
        !           249: end$

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