[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.3

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

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