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

Annotation of OpenXM/src/asir-contrib/testing/noro/ndbf.rr, Revision 1.4

1.1       noro        1: /* requires 'primdec' */
                      2:
                      3: #define TMP_H hhhhhhhh
                      4: #define TMP_S ssssssss
                      5: #define TMP_DS dssssssss
                      6: #define TMP_T t
                      7: #define TMP_DT dt
                      8: #define TMP_Y1 yyyyyyyy1
                      9: #define TMP_DY1 dyyyyyyyy1
                     10: #define TMP_Y2 yyyyyyyy2
                     11: #define TMP_DY2 dyyyyyyyy2
                     12:
                     13: if (!module_definedp("gr")) load("gr")$ else{ }$
                     14: if (!module_definedp("primdec")) load("primdec")$ else{ }$
                     15:   /* Empty for now. It will be used in a future. */
                     16:
                     17: /* toplevel */
                     18:
                     19: module ndbf$
                     20:
                     21: /* bfunction */
                     22:
                     23: localf bfunction, in_ww, in_ww_main, ann, ann_n$
                     24: localf ann0, psi, ww_weight, compare_first, generic_bfct$
                     25: localf generic_bfct_1, initial_part, bfct, indicial1, bfct_via_gbfct$
                     26: localf bfct_via_gbfct_weight, bfct_via_gbfct_weight_1, bfct_via_gbfct_weight_2$
                     27: localf weyl_minipolym, weyl_minipoly, weyl_nf, weyl_nf_quo_check$
                     28: localf weyl_nf_quo, weyl_nf_mod, b_subst, v_factorial, w_tdeg$
                     29: localf replace_vars_f, replace_vars_v, replace_var$
                     30: localf action_on_gfs, action_on_gfs_1$
1.2       noro       31: localf nd_gb_candidate$
1.1       noro       32:
                     33: /* stratification */
                     34:
                     35: localf weyl_subst, bfactor, gen_a, gen_d$
                     36: localf gen_dm, elimination, weyl_ideal_quotient, psi0$
                     37: localf bf_strat, bf_strat_stage2, bf_strat_stage3, bf_local$
                     38: localf conv_tdt, merge_tower, stratify_bf, tdt_homogenize$
                     39: localf sing, tower_in_p, subst_vars, compute_exponent$
                     40: localf zero_inclusion, weyl_divide_by_right, elim_mat, int_cons$
                     41: localf ideal_intersection$
                     42:
                     43: def bfunction(F)
                     44: {
                     45:        if ( type(Heu=getopt(heuristic)) == -1 ) Heu = 0;
                     46:        if ( type(Vord=getopt(vord)) == -1 || type(Vord) != 4 ) Vord = 0;
                     47:        if ( type(Wt=getopt(weight)) == -1 ) Wt = 0;
                     48:        L = in_ww(F|weight=Wt,heuristic=Heu,vord=Vord);
                     49:        Indata = L[0]; AllData = L[1]; VData = L[2];
                     50:        GIN = Indata[0]; VDV = Indata[1]; WVDV = AllData[4];
                     51:        W = Indata[4];
                     52:        dp_set_weight(W);
                     53:        B = weyl_minipoly(GIN,VDV,0,WVDV);
                     54:        dp_set_weight(0);
                     55:        return subst(B,s,-s-1);
                     56: }
                     57:
                     58: /*
                     59:        returns [InData,AllData,VData]
                     60:        InData = [GIN,VDV,V,DV,WtV]
                     61:        AllData = [G0,GIN0,VDV0,W,WVDV,WtV0]
                     62:        VData = [V0,DV0,T,DT]
                     63:        GIN0 = ini_(-W,W)(G0)
                     64:        WVDV = W[0]*V[0]*DV[0]+...
                     65: */
                     66:
                     67: def in_ww(F)
                     68: {
1.2       noro       69:        F = ptozp(F);
1.1       noro       70:        V = vars(F);
                     71:        N = length(V);
                     72:        D = newvect(N);
                     73:        Wt = getopt(weight);
                     74:        Vord = getopt(vord);
                     75:        if ( type(Wt) != 4 ) {
                     76:                if ( type(Vord) != 4 ) {
                     77:                for ( I = 0; I < N; I++ )
                     78:                        D[I] = [deg(F,V[I]),V[I]];
                     79:                qsort(D,compare_first);
                     80:                for ( V = [], I = 0; I < N; I++ )
                     81:                        V = cons(D[I][1],V);
                     82:                        V = reverse(V);
                     83:                } else
                     84:                        V = Vord;
                     85:                for ( I = 0, Wt = []; I < N; I++ )
                     86:                        Wt = cons(1,Wt);
                     87:        } else {
                     88:                Wt1 = vector(N);
                     89:                if ( type(Vord) != 4 ) {
                     90:                        for ( I = 0, F1 =F; I < N; I++ ) {
                     91:                                VI = Wt[2*I]; WI = Wt[2*I+1];
                     92:                                for ( J = 0; J < N; J++ )
                     93:                                        if ( VI == V[J] ) break;
                     94:                                F1 = subst(F1,VI,VI^WI);
                     95:                        }
                     96:                for ( I = 0; I < N; I++ )
                     97:                        D[I] = [deg(F1,V[I]),V[I]];
                     98:                qsort(D,compare_first);
                     99:                for ( V = [], I = 0; I < N; I++ )
                    100:                        V = cons(D[I][1],V);
                    101:                        V = reverse(V);
                    102:                } else
                    103:                        V = Vord;
                    104:                for ( I = 0; I < N; I++ ) {
                    105:                        VI = Wt[2*I]; WI = Wt[2*I+1];
                    106:                        for ( J = 0; J < N; J++ )
                    107:                                if ( VI == V[J] ) break;
                    108:                        Wt1[J] = WI;
                    109:                }
                    110:                Wt = vtol(Wt1);
                    111:        }
                    112:        Tdeg = w_tdeg(F,V,Wt);
                    113:        /* weight for [t,x1,...,xn,dt,dx1,...,dxn,h] */
                    114:        WtV = newvect(2*(N+1)+1);
                    115:        WtV[0] = Tdeg; WtV[N+1] = 1; WtV[2*(N+1)] = 1;
                    116:        /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */
                    117:        for ( I = 1; I <= N; I++ ) {
                    118:                WtV[I] = Wt[I-1]; WtV[N+1+I] = Tdeg-Wt[I-1]+1;
                    119:        }
                    120:        for ( I = N-1, DV = []; I >= 0; I-- )
                    121:                DV = cons(strtov("d"+rtostr(V[I])),DV);
                    122:
                    123:        B = [TMP_T-F];
                    124:        for ( I = 0; I < N; I++ ) {
                    125:                B = cons(DV[I]+diff(F,V[I])*TMP_DT,B);
                    126:        }
                    127:        V1 = cons(TMP_T,V); DV1 = cons(TMP_DT,DV);
                    128:        W = newvect(N+1); W[0] = 1;
                    129:     VW1 = [V1,DV1,WtV,W];
                    130:
                    131:        /* weight for [x1,...,xn,t,dx1,...,dxn,dt,h] */
                    132:        WtV = newvect(2*(N+1)+1);
                    133:        WtV[N] = Tdeg; WtV[2*N+1] = 1; WtV[2*(N+1)] = 1;
                    134:        for ( I = 0; I < N; I++ ) {
                    135:                WtV[I] = Wt[I]; WtV[N+1+I] = Tdeg-Wt[I]+1;
                    136:        }
                    137:        for ( I = N-1, DV = []; I >= 0; I-- )
                    138:                DV = cons(strtov("d"+rtostr(V[I])),DV);
                    139:
                    140:        B = [TMP_T-F];
                    141:        for ( I = 0; I < N; I++ ) {
                    142:                B = cons(DV[I]+diff(F,V[I])*TMP_DT,B);
                    143:        }
                    144:        V1 = append(V,[TMP_T]); DV1 = append(DV,[TMP_DT]);
                    145:        W = newvect(N+1); W[N] = 1;
                    146:     VW2 = [V1,DV1,WtV,W];
                    147:
                    148:        Heu = getopt(heuristic);
                    149:        if ( type(Heu) != -1 && Heu )
                    150:                L = in_ww_main(B,VW1,VW2);
                    151:        else
                    152:                L = in_ww_main(B,VW1,0);
                    153:        return append(L,[[V,DV,TMP_T,TMP_DT,Wt]]);
                    154: }
                    155:
                    156: /*
                    157:        returns [InData,AllData]
                    158:        InData = [GIN,VDV,V,DV,WtV]
                    159:        AllData = [G0,GIN0,VDV0,W,WVDV,WtV0]
                    160:        GIN0 = ini_(-W,W)(G0)
                    161:        WVDV = W[0]*V[0]*DV[0]+...
                    162: */
                    163:
                    164: def in_ww_main(F,VW1,VW2)
                    165: {
                    166:        V = VW1[0]; DV = VW1[1]; WtV = VW1[2]; W = VW1[3];
                    167:        dp_set_weight(WtV);
                    168:
                    169:        N = length(V);
                    170:        N2 = N*2;
                    171:
                    172:        /* If W is a list, convert it to a vector */
                    173:        if ( type(W) == 4 )
                    174:                W = newvect(length(W),W);
                    175:        dp_weyl_set_weight(W);
                    176:
                    177:        /* create a term order M in D<x,d> (DRL) */
                    178:        M = newmat(N2,N2);
                    179:        for ( J = 0; J < N2; J++ )
                    180:                M[0][J] = 1;
                    181:        for ( I = 1; I < N2; I++ )
                    182:                M[I][N2-I] = -1;
                    183:
                    184:        VDV = append(V,DV);
                    185:
                    186:        /* create a non-term order MW in D<x,d> */
                    187:        MW = newmat(N2+1,N2);
                    188:        for ( J = 0; J < N; J++ )
                    189:                MW[0][J] = -W[J];
                    190:        for ( ; J < N2; J++ )
                    191:                MW[0][J] = W[J-N];
                    192:        for ( I = 1; I <= N2; I++ )
                    193:                for ( J = 0; J < N2; J++ )
                    194:                        MW[I][J] = M[I-1][J];
                    195:
                    196:        /* create a homogenized term order MWH in D<x,d,h> */
                    197:        MWH = newmat(N2+2,N2+1);
                    198:        for ( J = 0; J <= N2; J++ )
                    199:                MWH[0][J] = 1;
                    200:        for ( I = 1; I <= N2+1; I++ )
                    201:                for ( J = 0; J < N2; J++ )
                    202:                        MWH[I][J] = MW[I-1][J];
                    203:
                    204:        /* homogenize F */
                    205:        VDVH = append(VDV,[TMP_H]);
                    206:        FH = map(dp_dtop,map(dp_homo,map(dp_ptod,F,VDV)),VDVH);
                    207:
1.2       noro      208: /*
                    209:  * FH is a GB w.r.t. any term order s.t. LT(FH)=[t,dx1,...,dxn]
                    210:  * Compute a groebner basis of FH w.r.t. MWH by modular change of
                    211:  * ordering.
                    212:  * Since F is Z-coef, LC(FH)=[1,...,1] and we can use any prime p
                    213:  * for trace algorithm.
                    214:  */
1.1       noro      215: /*     dp_gr_flags(["Top",1,"NoRA",1]); */
1.2       noro      216:        for ( I = 0, HC=[]; I <= N; I++ ) HC = cons(1,HC);
                    217:        GH = nd_gb_candidate(FH,VDVH,11,0,HC,1);
1.1       noro      218: /*     dp_gr_flags(["Top",0,"NoRA",0]); */
                    219:
                    220:        /* dehomigenize GH */
                    221:        G = map(subst,GH,TMP_H,1);
                    222:
                    223:        /* G is a groebner basis w.r.t. a non term order MW */
                    224:        /* take the initial part w.r.t. (-W,W) */
                    225:        GIN = map(initial_part,G,VDV,MW,W);
                    226:
                    227:        /* GIN is a groebner basis w.r.t. a term order M */
                    228:        /* As -W+W=0, gr_(-W,W)(D<x,d>) = D<x,d> */
                    229:
                    230:        /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */
                    231:        for ( I = 0, T = 0; I < N; I++ )
                    232:                T += W[I]*V[I]*DV[I];
                    233:
                    234:        AllData = [G,GIN,VDV,W,T,WtV];
                    235:        if ( VW2 ) {
1.2       noro      236:                /* take LC(GIN) w.r.t. DRL */
                    237:         dp_set_weight(WtV); dp_ord(0);
                    238:         HC = map(dp_hc,map(dp_ptod,GIN,VDV));
1.1       noro      239:                V2 = VW2[0]; DV2 = VW2[1]; WtV2 = VW2[2];
                    240:                VDV2 = append(V2,DV2);
                    241:                dp_set_weight(WtV2);
1.2       noro      242:                GIN2 = nd_gb_candidate(GIN,VDV2,0,0,HC,1);
1.1       noro      243:                InData = [GIN2,VDV2,V2,DV2,WtV2];
                    244:        } else {
                    245:                if ( 0 ) {
                    246:                        dp_set_weight(WtV);
                    247:                        GIN1 = nd_weyl_gr_postproc(GIN,VDV,0,0,0);
                    248:                        InData = [GIN1,VDV,V,DV,WtV];
                    249:                } else
                    250:                        InData = [GIN,VDV,V,DV,WtV];
                    251:        }
                    252:
                    253: /*     B = weyl_minipoly(GIN2,VDV2,0,T); */ /* M represents DRL order */
                    254:        WtV = dp_set_weight();
                    255:        dp_set_weight(0);
                    256:
                    257:        return [InData,AllData];
                    258: }
                    259:
                    260: /* annihilating ideal of F^s */
                    261:
                    262: def ann(F)
                    263: {
                    264:        if ( member(s,vars(F)) )
                    265:                error("ann : the variable 's' is reserved.");
1.4     ! noro      266:        F = ptozp(F);
1.1       noro      267:        V = vars(F);
                    268:        N = length(V);
                    269:        D = newvect(N);
1.4     ! noro      270:        if ( type(Wt=getopt(weight)) == -1 )
        !           271:                for ( I = N-1, Wt = []; I >= 0; I-- ) Wt = append([V[I],1],Wt);
1.1       noro      272:
1.4     ! noro      273:        Wt1 = vector(N);
        !           274:        for ( I = 0, F1 =F; I < N; I++ ) {
        !           275:                VI = Wt[2*I]; WI = Wt[2*I+1];
        !           276:                for ( J = 0; J < N; J++ )
        !           277:                        if ( VI == V[J] ) break;
        !           278:                F1 = subst(F1,VI,VI^WI);
        !           279:        }
        !           280:        for ( I = 0; I < N; I++ ) D[I] = [deg(F1,V[I]),V[I]];
1.1       noro      281:        qsort(D,compare_first);
1.4     ! noro      282:        for ( V = [], I = 0; I < N; I++ ) V = cons(D[I][1],V);
        !           283:        V = reverse(V);
        !           284:        for ( I = 0; I < N; I++ ) {
        !           285:                VI = Wt[2*I]; WI = Wt[2*I+1];
        !           286:                for ( J = 0; J < N; J++ ) if ( VI == V[J] ) break;
        !           287:                Wt1[J] = WI;
        !           288:        }
        !           289:        Wt = vtol(Wt1);
1.1       noro      290:
                    291:        for ( I = N-1, DV = []; I >= 0; I-- )
                    292:                DV = cons(strtov("d"+rtostr(V[I])),DV);
                    293:
                    294:        W = append([TMP_Y1,TMP_Y2,TMP_T],V);
                    295:        DW = append([TMP_DY1,TMP_DY2,TMP_DT],DV);
                    296:
                    297:        B = [1-TMP_Y1*TMP_Y2,TMP_T-TMP_Y1*F];
                    298:        for ( I = 0; I < N; I++ ) {
                    299:                B = cons(DV[I]+TMP_Y1*diff(F,V[I])*TMP_DT,B);
                    300:        }
                    301:
1.4     ! noro      302:        Tdeg = w_tdeg(F,V,Wt);
        !           303:        /* y1*y2-1, t-y1*f, dx1+y1*df/dx1*dt ... */
        !           304:        /* weight for [y1,y2,t,   x1,...,xn,  dy1,dy2, dt,dx1,...,dxn,   h]   */
        !           305:        /*              0  1 2    3      N3-1 N3  N3+1 N3+2              2*N3 */
        !           306:        /*              1  1 D+1  1      1    1   1    1  D       D      1    */
        !           307:        N3 = N+3;
        !           308:        WtV = newvect(2*N3+1);
        !           309:        WtV[0] = WtV[1] = 1; WtV[2] = Tdeg+1;
        !           310:        for ( I = 3; I <= N3+2; I++ ) WtV[I] = 1;
        !           311:        for ( ; I < 2*N3; I++ ) WtV[I] = Tdeg;
        !           312:        WtV[2*N3] = 1;
        !           313:
        !           314:        /* B is already a GB => modular change of ordering can be applied */
        !           315:        /* any prime is available => HC=[1] */
        !           316:        dp_set_weight(WtV);
        !           317:        G0 = nd_gb_candidate(B,append(W,DW),[[0,2],[0,length(W)*2-2]],0,[1],1);
        !           318:        dp_set_weight(0);
1.1       noro      319:        G1 = [];
                    320:        for ( T = G0; T != []; T = cdr(T) ) {
                    321:                E = car(T); VL = vars(E);
                    322:                if ( !member(TMP_Y1,VL) && !member(TMP_Y2,VL) )
                    323:                        G1 = cons(E,G1);
                    324:        }
                    325:        G2 = map(psi,G1,TMP_T,TMP_DT);
                    326:        G3 = map(subst,G2,TMP_T,-1-s);
                    327:        return G3;
                    328: }
                    329:
                    330: /* F = [F0,F1,...] */
                    331:
                    332: def ann_n(F)
                    333: {
                    334:        L = length(F);
                    335:        V = vars(F);
                    336:        N = length(V);
                    337:        D = newvect(N);
                    338:
                    339:        for ( I = N-1, DV = []; I >= 0; I-- )
                    340:                DV = cons(strtov("d"+rtostr(V[I])),DV);
                    341:        W = []; DW = [];
                    342:        for ( I = L-1; I >= 0; I-- ) {
                    343:                SI = rtostr(I);
                    344:                W = cons(strtov("t"+SI),W);
                    345:                DW = cons(strtov("dt"+SI),DW);
                    346:        }
                    347:        U = []; DU = [];
                    348:        for ( I = L-1; I >= 0; I-- ) {
                    349:                SI = rtostr(I);
                    350:                U = append([strtov("u"+SI),strtov("v"+SI)],U);
                    351:                DU = append([strtov("du"+SI),strtov("dv"+SI)],DU);
                    352:        }
                    353:
                    354:        B = [];
                    355:        for ( I = 0; I < N; I++ ) {
                    356:                T = DV[I];
                    357:                for ( J = 0; J < L; J++ )
                    358:                        T += U[2*J]*diff(F[J],V[I])*DW[J];
                    359:                B = cons(T,B);
                    360:        }
                    361:        for ( I = 0; I < L; I++ )
                    362:                B = append([W[I]-U[2*I]*F[I],1-U[2*I]*U[2*I+1]],B);
                    363:
                    364:        VA = append(U,append(W,V));
                    365:        DVA = append(DU,append(DW,DV));
                    366:        VDV = append(VA,DVA);
1.4     ! noro      367: #if 0
1.1       noro      368:        G0 = nd_weyl_gr(B,VDV,0,[[0,2*L],[0,length(VDV)-2*L]]);
1.4     ! noro      369: #else
        !           370:        G0 = nd_gb_candidate(B,VDV,[[0,2*L],[0,length(VDV)-2*L]],0,[1],1);
        !           371: #endif
1.1       noro      372:        G1 = [];
                    373:        for ( T = G0; T != []; T = cdr(T) ) {
                    374:                E = car(T); VL = vars(E);
                    375:                for ( TV = U; TV != []; TV = cdr(TV) )
                    376:                        if ( member(car(TV),VL) ) break;
                    377:                if ( TV == [] )
                    378:                        G1 = cons(E,G1);
                    379:        }
                    380:        G2 = G1;
                    381:        for ( I = 0; I < L; I++ ) {
                    382:                G2 = map(psi,G2,W[I],DW[I]);
                    383:                G2 = map(subst,G2,W[I],-1-strtov("s"+rtostr(I)));
                    384:        }
                    385:        return G2;
                    386: }
                    387:
                    388: /*
                    389:  * compute J_f|s=r, where r = the minimal integral root of global b_f(s)
                    390:  * ann0(F) returns [MinRoot,Ideal]
                    391:  */
                    392:
                    393: def ann0(F)
                    394: {
                    395:        F = subst(F,s,TMP_S);
                    396:        Ann = ann(F);
                    397:        Bf = bfunction(F);
                    398:
                    399:        FList = cdr(fctr(Bf));
                    400:        for ( T = FList, Min = 0; T != []; T = cdr(T) ) {
                    401:                LF = car(car(T));
                    402:                Root = -coef(LF,0)/coef(LF,1);
                    403:                if ( dn(Root) == 1 && Root < Min )
                    404:                        Min = Root;
                    405:        }
                    406:        return [Min,map(subst,Ann,s,Min,TMP_S,s,TMP_DS,ds)];
                    407: }
                    408:
                    409: def psi0(F,T,DT)
                    410: {
                    411:        D = dp_ptod(F,[T,DT]);
                    412:        Wmax = ww_weight(D);
                    413:        D1 = dp_rest(D);
                    414:        for ( ; D1; D1 = dp_rest(D1) )
                    415:                if ( ww_weight(D1) > Wmax )
                    416:                        Wmax = ww_weight(D1);
                    417:        for ( D1 = D, Dmax = 0; D1; D1 = dp_rest(D1) )
                    418:                if ( ww_weight(D1) == Wmax )
                    419:                        Dmax += dp_hm(D1);
                    420:        if ( Wmax >= 0 )
                    421:                Dmax = dp_weyl_mul(<<Wmax,0>>,Dmax);
                    422:        else
                    423:                Dmax = dp_weyl_mul(<<0,-Wmax>>,Dmax);
                    424:        Rmax = dp_dtop(Dmax,[T,DT]);
                    425:        return Rmax;
                    426: }
                    427:
                    428: def psi(F,T,DT)
                    429: {
                    430:        Rmax = psi0(F,T,DT);
                    431:        R = b_subst(subst(Rmax,DT,1),T);
                    432:        return R;
                    433: }
                    434:
                    435: def ww_weight(D)
                    436: {
                    437:        V = dp_etov(D);
                    438:        return V[1]-V[0];
                    439: }
                    440:
                    441: def compare_first(A,B)
                    442: {
                    443:        A0 = car(A);
                    444:        B0 = car(B);
                    445:        if ( A0 > B0 )
                    446:                return 1;
                    447:        else if ( A0 < B0 )
                    448:                return -1;
                    449:        else
                    450:                return 0;
                    451: }
                    452:
                    453: /* generic b-function w.r.t. weight vector W */
                    454:
                    455: def generic_bfct(F,V,DV,W)
                    456: {
                    457:        N = length(V);
                    458:        N2 = N*2;
                    459:
                    460:        /* If W is a list, convert it to a vector */
                    461:        if ( type(W) == 4 )
                    462:                W = newvect(length(W),W);
                    463:        dp_weyl_set_weight(W);
                    464:
                    465:        /* create a term order M in D<x,d> (DRL) */
                    466:        M = newmat(N2,N2);
                    467:        for ( J = 0; J < N2; J++ )
                    468:                M[0][J] = 1;
                    469:        for ( I = 1; I < N2; I++ )
                    470:                M[I][N2-I] = -1;
                    471:
                    472:        VDV = append(V,DV);
                    473:
                    474:        /* create a non-term order MW in D<x,d> */
                    475:        MW = newmat(N2+1,N2);
                    476:        for ( J = 0; J < N; J++ )
                    477:                MW[0][J] = -W[J];
                    478:        for ( ; J < N2; J++ )
                    479:                MW[0][J] = W[J-N];
                    480:        for ( I = 1; I <= N2; I++ )
                    481:                for ( J = 0; J < N2; J++ )
                    482:                        MW[I][J] = M[I-1][J];
                    483:
                    484:        /* create a homogenized term order MWH in D<x,d,h> */
                    485:        MWH = newmat(N2+2,N2+1);
                    486:        for ( J = 0; J <= N2; J++ )
                    487:                MWH[0][J] = 1;
                    488:        for ( I = 1; I <= N2+1; I++ )
                    489:                for ( J = 0; J < N2; J++ )
                    490:                        MWH[I][J] = MW[I-1][J];
                    491:
                    492:        /* homogenize F */
                    493:        VDVH = append(VDV,[TMP_H]);
                    494:        FH = map(dp_dtop,map(dp_homo,map(dp_ptod,F,VDV)),VDVH);
                    495:
                    496:        /* compute a groebner basis of FH w.r.t. MWH */
                    497:        dp_gr_flags(["Top",1,"NoRA",1]);
                    498:        GH = dp_weyl_gr_main(FH,VDVH,0,1,11);
                    499:        dp_gr_flags(["Top",0,"NoRA",0]);
                    500:
                    501:        /* dehomigenize GH */
                    502:        G = map(subst,GH,TMP_H,1);
                    503:
                    504:        /* G is a groebner basis w.r.t. a non term order MW */
                    505:        /* take the initial part w.r.t. (-W,W) */
                    506:        GIN = map(initial_part,G,VDV,MW,W);
                    507:
                    508:        /* GIN is a groebner basis w.r.t. a term order M */
                    509:        /* As -W+W=0, gr_(-W,W)(D<x,d>) = D<x,d> */
                    510:
                    511:        /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */
                    512:        for ( I = 0, T = 0; I < N; I++ )
                    513:                T += W[I]*V[I]*DV[I];
                    514:        B = weyl_minipoly(GIN,VDV,0,T); /* M represents DRL order */
                    515:        return B;
                    516: }
                    517:
                    518: /* all term reduction + interreduce */
                    519: def generic_bfct_1(F,V,DV,W)
                    520: {
                    521:        N = length(V);
                    522:        N2 = N*2;
                    523:
                    524:        /* If W is a list, convert it to a vector */
                    525:        if ( type(W) == 4 )
                    526:                W = newvect(length(W),W);
                    527:        dp_weyl_set_weight(W);
                    528:
                    529:        /* create a term order M in D<x,d> (DRL) */
                    530:        M = newmat(N2,N2);
                    531:        for ( J = 0; J < N2; J++ )
                    532:                M[0][J] = 1;
                    533:        for ( I = 1; I < N2; I++ )
                    534:                M[I][N2-I] = -1;
                    535:
                    536:        VDV = append(V,DV);
                    537:
                    538:        /* create a non-term order MW in D<x,d> */
                    539:        MW = newmat(N2+1,N2);
                    540:        for ( J = 0; J < N; J++ )
                    541:                MW[0][J] = -W[J];
                    542:        for ( ; J < N2; J++ )
                    543:                MW[0][J] = W[J-N];
                    544:        for ( I = 1; I <= N2; I++ )
                    545:                for ( J = 0; J < N2; J++ )
                    546:                        MW[I][J] = M[I-1][J];
                    547:
                    548:        /* create a homogenized term order MWH in D<x,d,h> */
                    549:        MWH = newmat(N2+2,N2+1);
                    550:        for ( J = 0; J <= N2; J++ )
                    551:                MWH[0][J] = 1;
                    552:        for ( I = 1; I <= N2+1; I++ )
                    553:                for ( J = 0; J < N2; J++ )
                    554:                        MWH[I][J] = MW[I-1][J];
                    555:
                    556:        /* homogenize F */
                    557:        VDVH = append(VDV,[TMP_H]);
                    558:        FH = map(dp_dtop,map(dp_homo,map(dp_ptod,F,VDV)),VDVH);
                    559:
                    560:        /* compute a groebner basis of FH w.r.t. MWH */
                    561: /*     dp_gr_flags(["Top",1,"NoRA",1]); */
                    562:        GH = dp_weyl_gr_main(FH,VDVH,0,1,11);
                    563: /*     dp_gr_flags(["Top",0,"NoRA",0]); */
                    564:
                    565:        /* dehomigenize GH */
                    566:        G = map(subst,GH,TMP_H,1);
                    567:
                    568:        /* G is a groebner basis w.r.t. a non term order MW */
                    569:        /* take the initial part w.r.t. (-W,W) */
                    570:        GIN = map(initial_part,G,VDV,MW,W);
                    571:
                    572:        /* GIN is a groebner basis w.r.t. a term order M */
                    573:        /* As -W+W=0, gr_(-W,W)(D<x,d>) = D<x,d> */
                    574:
                    575:        /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */
                    576:        for ( I = 0, T = 0; I < N; I++ )
                    577:                T += W[I]*V[I]*DV[I];
                    578:        B = weyl_minipoly(GIN,VDV,0,T); /* M represents DRL order */
                    579:        return B;
                    580: }
                    581:
                    582: def initial_part(F,V,MW,W)
                    583: {
                    584:        N2 = length(V);
                    585:        N = N2/2;
                    586:        dp_ord(MW);
                    587:        DF = dp_ptod(F,V);
                    588:        R = dp_hm(DF);
                    589:        DF = dp_rest(DF);
                    590:
                    591:        E = dp_etov(R);
                    592:        for ( I = 0, TW = 0; I < N; I++ )
                    593:                TW += W[I]*(-E[I]+E[N+I]);
                    594:        RW = TW;
                    595:
                    596:        for ( ; DF; DF = dp_rest(DF) ) {
                    597:                E = dp_etov(DF);
                    598:                for ( I = 0, TW = 0; I < N; I++ )
                    599:                        TW += W[I]*(-E[I]+E[N+I]);
                    600:                if ( TW == RW )
                    601:                        R += dp_hm(DF);
                    602:                else if ( TW < RW )
                    603:                        break;
                    604:                else
                    605:                        error("initial_part : cannot happen");
                    606:        }
                    607:        return dp_dtop(R,V);
                    608:
                    609: }
                    610:
                    611: /* b-function of F ? */
                    612:
                    613: def bfct(F)
                    614: {
                    615:        /* XXX */
                    616:        F = replace_vars_f(F);
                    617:
                    618:        V = vars(F);
                    619:        N = length(V);
                    620:        D = newvect(N);
                    621:
                    622:        for ( I = 0; I < N; I++ )
                    623:                D[I] = [deg(F,V[I]),V[I]];
                    624:        qsort(D,compare_first);
                    625:        for ( V = [], I = 0; I < N; I++ )
                    626:                V = cons(D[I][1],V);
                    627:        for ( I = N-1, DV = []; I >= 0; I-- )
                    628:                DV = cons(strtov("d"+rtostr(V[I])),DV);
                    629:        V1 = cons(s,V); DV1 = cons(ds,DV);
                    630:
                    631:        G0 = indicial1(F,reverse(V));
                    632:        G1 = dp_weyl_gr_main(G0,append(V1,DV1),0,1,0);
                    633:        Minipoly = weyl_minipoly(G1,append(V1,DV1),0,s);
                    634:        return Minipoly;
                    635: }
                    636:
                    637: /* called from bfct() only */
                    638:
                    639: def indicial1(F,V)
                    640: {
                    641:        W = append([y1,t],V);
                    642:        N = length(V);
                    643:        B = [t-y1*F];
                    644:        for ( I = N-1, DV = []; I >= 0; I-- )
                    645:                DV = cons(strtov("d"+rtostr(V[I])),DV);
                    646:        DW = append([dy1,dt],DV);
                    647:        for ( I = 0; I < N; I++ ) {
                    648:                B = cons(DV[I]+y1*diff(F,V[I])*dt,B);
                    649:        }
                    650:        dp_nelim(1);
                    651:
                    652:        /* homogenized (heuristics) */
                    653:        G0 = dp_weyl_gr_main(B,append(W,DW),1,0,6);
                    654:        G1 = map(subst,G0,y1,1);
                    655:        G2 = map(psi,G1,t,dt);
                    656:        G3 = map(subst,G2,t,-s-1);
                    657:        return G3;
                    658: }
                    659:
                    660: /* b-function computation via generic_bfct() (experimental) */
                    661:
                    662: def bfct_via_gbfct(F)
                    663: {
                    664:        V = vars(F);
                    665:        N = length(V);
                    666:        D = newvect(N);
                    667:
                    668:        for ( I = 0; I < N; I++ )
                    669:                D[I] = [deg(F,V[I]),V[I]];
                    670:        qsort(D,compare_first);
                    671:        for ( V = [], I = 0; I < N; I++ )
                    672:                V = cons(D[I][1],V);
                    673:        V = reverse(V);
                    674:        for ( I = N-1, DV = []; I >= 0; I-- )
                    675:                DV = cons(strtov("d"+rtostr(V[I])),DV);
                    676:
                    677:        B = [TMP_T-F];
                    678:        for ( I = 0; I < N; I++ ) {
                    679:                B = cons(DV[I]+diff(F,V[I])*TMP_DT,B);
                    680:        }
                    681:        V1 = cons(TMP_T,V); DV1 = cons(TMP_DT,DV);
                    682:        W = newvect(N+1);
                    683:        W[0] = 1;
                    684:        R = generic_bfct(B,V1,DV1,W);
                    685:
                    686:        return subst(R,s,-s-1);
                    687: }
                    688:
                    689: /* use an order s.t. [t,x,y,z,...,dt,dx,dy,dz,...,h] */
                    690:
                    691: def bfct_via_gbfct_weight(F,V)
                    692: {
                    693:        N = length(V);
                    694:        D = newvect(N);
                    695:        Wt = getopt(weight);
                    696:        if ( type(Wt) != 4 ) {
                    697:                for ( I = 0, Wt = []; I < N; I++ )
                    698:                        Wt = cons(1,Wt);
                    699:        }
                    700:        Tdeg = w_tdeg(F,V,Wt);
                    701:        WtV = newvect(2*(N+1)+1);
                    702:        WtV[0] = Tdeg;
                    703:        WtV[N+1] = 1;
                    704:        /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */
                    705:        for ( I = 1; I <= N; I++ ) {
                    706:                WtV[I] = Wt[I-1];
                    707:                WtV[N+1+I] = Tdeg-Wt[I-1]+1;
                    708:        }
                    709:        WtV[2*(N+1)] = 1;
                    710:        dp_set_weight(WtV);
                    711:        for ( I = N-1, DV = []; I >= 0; I-- )
                    712:                DV = cons(strtov("d"+rtostr(V[I])),DV);
                    713:
                    714:        B = [TMP_T-F];
                    715:        for ( I = 0; I < N; I++ ) {
                    716:                B = cons(DV[I]+diff(F,V[I])*TMP_DT,B);
                    717:        }
                    718:        V1 = cons(TMP_T,V); DV1 = cons(TMP_DT,DV);
                    719:        W = newvect(N+1);
                    720:        W[0] = 1;
                    721:        R = generic_bfct_1(B,V1,DV1,W);
                    722:        dp_set_weight(0);
                    723:        return subst(R,s,-s-1);
                    724: }
                    725:
                    726: /* use an order s.t. [x,y,z,...,t,dx,dy,dz,...,dt,h] */
                    727:
                    728: def bfct_via_gbfct_weight_1(F,V)
                    729: {
                    730:        N = length(V);
                    731:        D = newvect(N);
                    732:        Wt = getopt(weight);
                    733:        if ( type(Wt) != 4 ) {
                    734:                for ( I = 0, Wt = []; I < N; I++ )
                    735:                        Wt = cons(1,Wt);
                    736:        }
                    737:        Tdeg = w_tdeg(F,V,Wt);
                    738:        WtV = newvect(2*(N+1)+1);
                    739:        /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */
                    740:        for ( I = 0; I < N; I++ ) {
                    741:                WtV[I] = Wt[I];
                    742:                WtV[N+1+I] = Tdeg-Wt[I]+1;
                    743:        }
                    744:        WtV[N] = Tdeg;
                    745:        WtV[2*N+1] = 1;
                    746:        WtV[2*(N+1)] = 1;
                    747:        dp_set_weight(WtV);
                    748:        for ( I = N-1, DV = []; I >= 0; I-- )
                    749:                DV = cons(strtov("d"+rtostr(V[I])),DV);
                    750:
                    751:        B = [TMP_T-F];
                    752:        for ( I = 0; I < N; I++ ) {
                    753:                B = cons(DV[I]+diff(F,V[I])*TMP_DT,B);
                    754:        }
                    755:        V1 = append(V,[TMP_T]); DV1 = append(DV,[TMP_DT]);
                    756:        W = newvect(N+1);
                    757:        W[N] = 1;
                    758:        R = generic_bfct_1(B,V1,DV1,W);
                    759:        dp_set_weight(0);
                    760:        return subst(R,s,-s-1);
                    761: }
                    762:
                    763: def bfct_via_gbfct_weight_2(F,V)
                    764: {
                    765:        N = length(V);
                    766:        D = newvect(N);
                    767:        Wt = getopt(weight);
                    768:        if ( type(Wt) != 4 ) {
                    769:                for ( I = 0, Wt = []; I < N; I++ )
                    770:                        Wt = cons(1,Wt);
                    771:        }
                    772:        Tdeg = w_tdeg(F,V,Wt);
                    773:
                    774:        /* a weight for the first GB computation */
                    775:        /* [t,x1,...,xn,dt,dx1,...,dxn,h] */
                    776:        WtV = newvect(2*(N+1)+1);
                    777:        WtV[0] = Tdeg;
                    778:        WtV[N+1] = 1;
                    779:        WtV[2*(N+1)] = 1;
                    780:        /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */
                    781:        for ( I = 1; I <= N; I++ ) {
                    782:                WtV[I] = Wt[I-1];
                    783:                WtV[N+1+I] = Tdeg-Wt[I-1]+1;
                    784:        }
                    785:        dp_set_weight(WtV);
                    786:
                    787:        /* a weight for the second GB computation */
                    788:        /* [x1,...,xn,t,dx1,...,dxn,dt,h] */
                    789:        WtV2 = newvect(2*(N+1)+1);
                    790:        WtV2[N] = Tdeg;
                    791:        WtV2[2*N+1] = 1;
                    792:        WtV2[2*(N+1)] = 1;
                    793:        for ( I = 0; I < N; I++ ) {
                    794:                WtV2[I] = Wt[I];
                    795:                WtV2[N+1+I] = Tdeg-Wt[I]+1;
                    796:        }
                    797:
                    798:        for ( I = N-1, DV = []; I >= 0; I-- )
                    799:                DV = cons(strtov("d"+rtostr(V[I])),DV);
                    800:
                    801:        B = [TMP_T-F];
                    802:        for ( I = 0; I < N; I++ ) {
                    803:                B = cons(DV[I]+diff(F,V[I])*TMP_DT,B);
                    804:        }
                    805:        V1 = cons(TMP_T,V); DV1 = cons(TMP_DT,DV);
                    806:        V2 = append(V,[TMP_T]); DV2 = append(DV,[TMP_DT]);
                    807:        W = newvect(N+1,[1]);
                    808:        dp_weyl_set_weight(W);
                    809:
                    810:        VDV = append(V1,DV1);
                    811:        N1 = length(V1);
                    812:        N2 = N1*2;
                    813:
                    814:        /* create a non-term order MW in D<x,d> */
                    815:        MW = newmat(N2+1,N2);
                    816:        for ( J = 0; J < N1; J++ ) {
                    817:                MW[0][J] = -W[J]; MW[0][N1+J] = W[J];
                    818:        }
                    819:        for ( J = 0; J < N2; J++ ) MW[1][J] = 1;
                    820:        for ( I = 2; I <= N2; I++ ) MW[I][N2-I+1] = -1;
                    821:
                    822:        /* homogenize F */
                    823:        VDVH = append(VDV,[TMP_H]);
                    824:        FH = map(dp_dtop,map(dp_homo,map(dp_ptod,B,VDV)),VDVH);
                    825:
                    826:        /* compute a groebner basis of FH w.r.t. MWH */
                    827:        GH = dp_weyl_gr_main(FH,VDVH,0,1,11);
                    828:
                    829:        /* dehomigenize GH */
                    830:        G = map(subst,GH,TMP_H,1);
                    831:
                    832:        /* G is a groebner basis w.r.t. a non term order MW */
                    833:        /* take the initial part w.r.t. (-W,W) */
                    834:        GIN = map(initial_part,G,VDV,MW,W);
                    835:
                    836:        /* GIN is a groebner basis w.r.t. a term order M */
                    837:        /* As -W+W=0, gr_(-W,W)(D<x,d>) = D<x,d> */
                    838:
                    839:        /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */
                    840:        for ( I = 0, T = 0; I < N1; I++ )
                    841:                T += W[I]*V1[I]*DV1[I];
                    842:
                    843:        /* change of ordering from VDV to VDV2 */
                    844:        VDV2 = append(V2,DV2);
                    845:        dp_set_weight(WtV2);
                    846:        for ( Pind = 0; ; Pind++ ) {
                    847:                Prime = lprime(Pind);
                    848:                GIN2 = dp_weyl_gr_main(GIN,VDV2,0,-Prime,0);
                    849:                if ( GIN2 ) break;
                    850:        }
                    851:
                    852:        R = weyl_minipoly(GIN2,VDV2,0,T); /* M represents DRL order */
                    853:        dp_set_weight(0);
                    854:        return subst(R,s,-s-1);
                    855: }
                    856:
                    857: /* minimal polynomial of s; modular computation */
                    858:
                    859: def weyl_minipolym(G,V,O,M,V0)
                    860: {
                    861:        N = length(V);
                    862:        Len = length(G);
                    863:        dp_ord(O);
                    864:        setmod(M);
                    865:        PS = newvect(Len);
                    866:        PS0 = newvect(Len);
                    867:
                    868:        for ( I = 0, T = G; T != []; T = cdr(T), I++ )
                    869:                PS0[I] = dp_ptod(car(T),V);
                    870:        for ( I = 0, T = G; T != []; T = cdr(T), I++ )
                    871:                PS[I] = dp_mod(dp_ptod(car(T),V),M,[]);
                    872:
                    873:        for ( I = Len - 1, GI = []; I >= 0; I-- )
                    874:                GI = cons(I,GI);
                    875:
                    876:        U = dp_mod(dp_ptod(V0,V),M,[]);
                    877:        U = dp_weyl_nf_mod(GI,U,PS,1,M);
                    878:
                    879:        T = dp_mod(<<0>>,M,[]);
                    880:        TT = dp_mod(dp_ptod(1,V),M,[]);
                    881:        G = H = [[TT,T]];
                    882:
                    883:        for ( I = 1; ; I++ ) {
                    884:                if ( dp_gr_print() )
                    885:                        print(".",2);
                    886:                T = dp_mod(<<I>>,M,[]);
                    887:
                    888:                TT = dp_weyl_nf_mod(GI,dp_weyl_mul_mod(TT,U,M),PS,1,M);
                    889:                H = cons([TT,T],H);
                    890:                L = dp_lnf_mod([TT,T],G,M);
                    891:                if ( !L[0] ) {
                    892:                        if ( dp_gr_print() )
                    893:                                print("");
                    894:                        return dp_dtop(L[1],[t]); /* XXX */
                    895:                } else
                    896:                        G = insert(G,L);
                    897:        }
                    898: }
                    899:
                    900: /* minimal polynomial of s over Q */
                    901:
                    902: def weyl_minipoly(G0,V0,O0,P)
                    903: {
                    904:        HM = hmlist(G0,V0,O0);
                    905:
                    906:        N = length(V0);
                    907:        Len = length(G0);
                    908:        dp_ord(O0);
                    909:        PS = newvect(Len);
                    910:        for ( I = 0, T = G0, HL = []; T != []; T = cdr(T), I++ )
                    911:                PS[I] = dp_ptod(car(T),V0);
                    912:        for ( I = Len - 1, GI = []; I >= 0; I-- )
                    913:                GI = cons(I,GI);
                    914:        PSM = newvect(Len);
                    915:        DP = dp_ptod(P,V0);
                    916:
                    917:        for ( Pind = 0; ; Pind++ ) {
                    918:                Prime = lprime(Pind);
                    919:                if ( !valid_modulus(HM,Prime) )
                    920:                        continue;
                    921:                setmod(Prime);
                    922:                for ( I = 0, T = G0, HL = []; T != []; T = cdr(T), I++ )
                    923:                        PSM[I] = dp_mod(dp_ptod(car(T),V0),Prime,[]);
                    924:
                    925:                NFP = weyl_nf(GI,DP,1,PS);
                    926:                NFPM = dp_mod(NFP[0],Prime,[])/ptomp(NFP[1],Prime);
                    927:
                    928:                NF = [[dp_ptod(1,V0),1]];
                    929:                LCM = 1;
                    930:
                    931:                TM = dp_mod(<<0>>,Prime,[]);
                    932:                TTM = dp_mod(dp_ptod(1,V0),Prime,[]);
                    933:                GM = NFM = [[TTM,TM]];
                    934:
                    935:                for ( D = 1; ; D++ ) {
                    936:                        if ( dp_gr_print() )
                    937:                                print(".",2);
                    938:                        NFPrev = car(NF);
                    939:                        NFJ = weyl_nf(GI,
                    940:                                dp_weyl_mul(NFP[0],NFPrev[0]),NFP[1]*NFPrev[1],PS);
                    941:                        NFJ = remove_cont(NFJ);
                    942:                        NF = cons(NFJ,NF);
                    943:                        LCM = ilcm(LCM,NFJ[1]);
                    944:
                    945:                        /* modular computation */
                    946:                        TM = dp_mod(<<D>>,Prime,[]);
                    947:                        TTM = dp_mod(NFJ[0],Prime,[])/ptomp(NFJ[1],Prime);
                    948:                        NFM = cons([TTM,TM],NFM);
                    949:                        LM = dp_lnf_mod([TTM,TM],GM,Prime);
                    950:                        if ( !LM[0] )
                    951:                                break;
                    952:                        else
                    953:                                GM = insert(GM,LM);
                    954:                }
                    955:
                    956:                if ( dp_gr_print() )
                    957:                        print("");
                    958:                U = NF[0][0]*idiv(LCM,NF[0][1]);
                    959:                Coef = [];
                    960:                for ( J = D-1; J >= 0; J-- ) {
                    961:                        Coef = cons(strtov("u"+rtostr(J)),Coef);
                    962:                        U += car(Coef)*NF[D-J][0]*idiv(LCM,NF[D-J][1]);
                    963:                }
                    964:
                    965:                for ( UU = U, Eq = []; UU; UU = dp_rest(UU) )
                    966:                        Eq = cons(dp_hc(UU),Eq);
                    967:                M = etom([Eq,Coef]);
                    968:                B = henleq(M,Prime);
                    969:                if ( dp_gr_print() )
                    970:                        print("");
                    971:                if ( B ) {
                    972:                        R = 0;
                    973:                        for ( I = 0; I < D; I++ )
                    974:                                R += B[0][I]*s^I;
                    975:                        R += B[1]*s^D;
                    976:                        return R;
                    977:                }
                    978:        }
                    979: }
                    980:
                    981: def weyl_nf(B,G,M,PS)
                    982: {
                    983:        for ( D = 0; G; ) {
                    984:                for ( U = 0, L = B; L != []; L = cdr(L) ) {
                    985:                        if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
                    986:                                GCD = igcd(dp_hc(G),dp_hc(R));
                    987:                                CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
                    988:                                U = CG*G-dp_weyl_mul(CR*dp_subd(G,R),R);
                    989:                                if ( !U )
                    990:                                        return [D,M];
                    991:                                D *= CG; M *= CG;
                    992:                                break;
                    993:                        }
                    994:                }
                    995:                if ( U )
                    996:                        G = U;
                    997:                else {
                    998:                        D += dp_hm(G); G = dp_rest(G);
                    999:                }
                   1000:        }
                   1001:        return [D,M];
                   1002: }
                   1003:
                   1004: def weyl_nf_quo_check(G,PS,R)
                   1005: {
                   1006:        D = R[0]; M = R[1]; Coef = R[2];
                   1007:        Len = length(PS);
                   1008:        T = 0;
                   1009:        for ( I = 0; I < Len; I++ )
                   1010:                T += dp_weyl_mul(Coef[I],PS[I]);
                   1011:        return (M*G-T)==D;
                   1012: }
                   1013:
                   1014: def weyl_nf_quo(B,G,M,PS)
                   1015: {
                   1016:        Len = length(PS);
                   1017:        Coef = vector(Len);
                   1018:        for ( D = 0; G; ) {
                   1019:                for ( U = 0, L = B; L != []; L = cdr(L) ) {
                   1020:                        if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
                   1021:                                GCD = igcd(dp_hc(G),dp_hc(R));
                   1022:                                CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
                   1023:                                for ( I = 0; I < Len; I++ ) Coef[I] *= CG;
                   1024:                                Q = CR*dp_subd(G,R);
                   1025:                                Coef[car(L)] += Q;
                   1026:                                U = CG*G-dp_weyl_mul(Q,R);
                   1027:                                D *= CG; M *= CG;
                   1028:                                if ( !U )
                   1029:                                        return [D,M,Coef];
                   1030:                                break;
                   1031:                        }
                   1032:                }
                   1033:                if ( U )
                   1034:                        G = U;
                   1035:                else {
                   1036:                        D += dp_hm(G); G = dp_rest(G);
                   1037:                }
                   1038:        }
                   1039:        return [D,M,Coef];
                   1040: }
                   1041:
                   1042: def weyl_nf_mod(B,G,PS,Mod)
                   1043: {
                   1044:        for ( D = 0; G; ) {
                   1045:                for ( U = 0, L = B; L != []; L = cdr(L) ) {
                   1046:                        if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
                   1047:                                CR = dp_hc(G)/dp_hc(R);
                   1048:                                U = G-dp_weyl_mul_mod(CR*dp_mod(dp_subd(G,R),Mod,[]),R,Mod);
                   1049:                                if ( !U )
                   1050:                                        return D;
                   1051:                                break;
                   1052:                        }
                   1053:                }
                   1054:                if ( U )
                   1055:                        G = U;
                   1056:                else {
                   1057:                        D += dp_hm(G); G = dp_rest(G);
                   1058:                }
                   1059:        }
                   1060:        return D;
                   1061: }
                   1062:
                   1063: def b_subst(F,V)
                   1064: {
                   1065:        D = deg(F,V);
                   1066:        C = newvect(D+1);
                   1067:        for ( I = D; I >= 0; I-- )
                   1068:                C[I] = coef(F,I,V);
                   1069:        for ( I = 0, R = 0; I <= D; I++ )
                   1070:                if ( C[I] )
                   1071:                        R += C[I]*v_factorial(V,I);
                   1072:        return R;
                   1073: }
                   1074:
                   1075: def v_factorial(V,N)
                   1076: {
                   1077:        for ( J = N-1, R = 1; J >= 0; J-- )
                   1078:                R *= V-J;
                   1079:        return R;
                   1080: }
                   1081:
                   1082: def w_tdeg(F,V,W)
                   1083: {
                   1084:        dp_set_weight(newvect(length(W),W));
                   1085:        T = dp_ptod(F,V);
                   1086:        for ( R = 0; T; T = cdr(T) ) {
                   1087:                D = dp_td(T);
                   1088:                if ( D > R ) R = D;
                   1089:        }
                   1090:        return R;
                   1091: }
                   1092:
                   1093: def replace_vars_f(F)
                   1094: {
                   1095:        return subst(F,s,TMP_S,t,TMP_T,y1,TMP_Y1,y2,TMP_Y2);
                   1096: }
                   1097:
                   1098: def replace_vars_v(V)
                   1099: {
                   1100:        V = replace_var(V,s,TMP_S);
                   1101:        V = replace_var(V,t,TMP_T);
                   1102:        V = replace_var(V,y1,TMP_Y1);
                   1103:        V = replace_var(V,y2,TMP_Y2);
                   1104:        return V;
                   1105: }
                   1106:
                   1107: def replace_var(V,X,Y)
                   1108: {
                   1109:        V = reverse(V);
                   1110:        for ( R = []; V != []; V = cdr(V) ) {
                   1111:                Z = car(V);
                   1112:                if ( Z == X )
                   1113:                        R = cons(Y,R);
                   1114:                else
                   1115:                        R = cons(Z,R);
                   1116:        }
                   1117:        return R;
                   1118: }
                   1119:
                   1120: def action_on_gfs(P,V,GFS)
                   1121: {
                   1122:        DP = dp_ptod(P,V);
                   1123:        N = length(V)/2;
                   1124:        for ( I = N-1, V0 = []; I >= 0; I-- )
                   1125:                V0 = cons(V[I],V0);
                   1126:        R = [];
                   1127:        for ( T = DP; T; T = dp_rest(T) )
                   1128:                R = cons(action_on_gfs_1(dp_hm(T),N,V0,GFS),R);
                   1129:        D = coef(car(R)[2],0);
                   1130:        for ( T = cdr(R); T != []; T = cdr(T) ) {
                   1131:                Di = coef(car(T)[2],0);
                   1132:                if ( Di < D )
                   1133:                        D = Di;
                   1134:        }
                   1135:        F = GFS[1];
                   1136:        for ( T = R, G = 0; T != []; T = cdr(T) )
                   1137:                G += car(T)[0]*F^(car(T)[2]-(s+D));
                   1138:        while ( 1 ) {
                   1139:                G1 = tdiv(G,F);
                   1140:                if ( G1 ) {
                   1141:                        G = G1;
                   1142:                        D++;
                   1143:                } else
                   1144:                        return [G,F,s+D];
                   1145:        }
                   1146: }
                   1147:
                   1148: def action_on_gfs_1(M,N,V,GFS)
                   1149: {
                   1150:        G = GFS[0];
                   1151:        F = GFS[1];
                   1152:        S = GFS[2];
                   1153:        C = dp_hc(M);
                   1154:        E = dp_etov(M);
                   1155:        for ( I = 0; I < N; I++ ) {
                   1156:                VI = V[I];
                   1157:                C *= VI^E[I];
                   1158:                DFVI = diff(F,VI);
                   1159:                for ( J = 0, EI = E[I+N]; J < EI; J++, S-- )
                   1160:                        G = diff(G,VI)*F+S*G*DFVI;
                   1161:        }
                   1162:        return [C*G,F,S];
                   1163: }
                   1164:
                   1165: /* stratification */
                   1166:
                   1167: def weyl_subst(F,P,V)
                   1168: {
                   1169:        VF = var(F);
                   1170:        D = deg(F,VF);
                   1171:        P = dp_ptod(P,V);
                   1172:        One = dp_ptod(1,V);
                   1173:        for ( R = 0, I = D; I >= 0; I-- )
                   1174:                R = dp_weyl_mul(R,P)+coef(F,I,VF)*One;
                   1175:        return dp_dtop(R,V);
                   1176: }
                   1177:
                   1178: def bfactor(F)
                   1179: {
                   1180:        L=length(F);
                   1181:        for(I=0,B=1;I<L;I++)B*=F[I][0]^F[I][1];
                   1182:        return fctr(B);
                   1183: }
                   1184:
                   1185: def gen_a(K)
                   1186: {
                   1187:        D = x^(K+1);
                   1188:        W = [];
                   1189:        for ( I = 1; I <= K; I++ ) {
                   1190:                D += (V=strtov("u"+rtostr(K-I+1)))*x^(K-I);
                   1191:                W = cons(I+1,cons(V,W));
                   1192:        }
                   1193:        F = res(x,D,diff(D,x));
                   1194:        return [D,F,reverse(W)];
                   1195: }
                   1196:
                   1197: def gen_d(K)
                   1198: {
                   1199:        D = x^2*y+y^(K-1)+u1+u2*x+u3*x^2;
                   1200:        W = reverse([u1,2*K-2,u2,K,u3,2]);
                   1201:        U = [u3,u2,u1];
                   1202:        for ( I = 4; I <= K; I++ ) {
                   1203:                D += (V=strtov("u"+rtostr(I)))*y^(I-3);
                   1204:                W = cons((2*K-2)-2*(I-3),cons(V,W));
                   1205:                U = cons(V,U);
                   1206:        }
                   1207:        B = [D,diff(D,x),diff(D,y)];
                   1208:        G = nd_gr_trace(B,append([x,y],U),1,1,0);
                   1209:        G = nd_gr_trace(G,append([x,y],U),1,-1,[[0,2],[0,K]]);
                   1210:        E = elimination(G,U);
                   1211:        F = E[0];
                   1212:        return [D,F,reverse(W)];
                   1213: }
                   1214:
                   1215: def gen_dm(K)
                   1216: {
                   1217:        D = x^2*y-y^(K-1)+u1+u2*x+u3*x^2;
                   1218:        W = reverse([u1,2*K-2,u2,K,u3,2]);
                   1219:        U = [u3,u2,u1];
                   1220:        for ( I = 4; I <= K; I++ ) {
                   1221:                D += (V=strtov("u"+rtostr(I)))*y^(I-3);
                   1222:                W = cons((2*K-2)-2*(I-3),cons(V,W));
                   1223:                U = cons(V,U);
                   1224:        }
                   1225:        B = [D,diff(D,x),diff(D,y)];
                   1226:        G = nd_gr_trace(B,append([x,y],U),1,1,0);
                   1227:        G = nd_gr_trace(G,append([x,y],U),1,-1,[[0,2],[0,K]]);
                   1228:        E = elimination(G,U);
                   1229:        F = E[0];
                   1230:        return [D,F,reverse(W)];
                   1231: }
                   1232:
                   1233: def elimination(G,V)
                   1234:        {
                   1235:        ANS=[];
                   1236:        NG=length(G);
                   1237:
                   1238:        for (I=NG-1;I>=0;I--)
                   1239:                {
                   1240:                VSet=vars(G[I]);
                   1241:                DIFF=setminus(VSet,V);
                   1242:
                   1243:                if ( DIFF ==[] )
                   1244:                        {
                   1245:                        ANS=cons(G[I],ANS);
                   1246:                        }
                   1247:                }
                   1248:        return ANS;
                   1249:        }
                   1250:
                   1251: def weyl_ideal_quotient(B,F,VDV)
                   1252: {
                   1253:        T = ttttt; DT = dttttt;
                   1254:        J = cons((1-T)*F,vtol(ltov(B)*T));
                   1255:        N = length(VDV); N1 = N/2;
                   1256:        for ( I = N1-1, V1 = []; I >= 0; I-- )
                   1257:                V1 = cons(VDV[I],V1);
                   1258:        for ( I = 0, VDV1 = VDV; I < N1; I++ ) VDV1 = cdr(VDV1);
                   1259:        VDV1 = append(cons(T,V1),cons(DT,VDV1));
                   1260:        O1 = [[0,1],[0,N+1]];
                   1261:        GJ = nd_weyl_gr(J,VDV1,0,O1);
                   1262:        R = elimination(GJ,VDV);
                   1263:        return R;
                   1264: }
                   1265:
                   1266: def bf_strat(F)
                   1267: {
                   1268:        dp_ord(0);
                   1269:        T0 = time();
                   1270:        if ( type(Heu=getopt(heuristic)) == -1 ) Heu = 0;
                   1271:        if ( type(Vord=getopt(vord)) == -1 || type(Vord) != 4 ) Vord = 0;
                   1272:        if ( type(Wt=getopt(weight)) == -1 ) Wt = 0;
                   1273:        L = in_ww(F|weight=Wt,heuristic=Heu,vord=Vord);
                   1274:        T1 = time();
                   1275:        print(["in_ww",(T1[0]+T1[1])-(T0[0]+T0[1])]);
                   1276:
                   1277:        /* shortcuts */
                   1278:        V = vars(F);
                   1279:        dp_set_weight(0);
                   1280:        dp_ord(0);
                   1281:        Sing = sing(F,V);
                   1282:        if ( Sing[0] == 1 || Sing[0] == -1 ) {
                   1283:                return [[[F],[1],[[s+1,1]]],[[0],[F],[]]];
                   1284:        } else if ( zero_dim(Sing,V,0) ) {
                   1285:                N = length(V);
                   1286:                P0 = [];
                   1287:                for ( I = 0; I < N; I++ ) {
                   1288:                        M = minipoly(Sing,V,0,V[I],TMP_S);
                   1289:                        MF = cdr(fctr(M));
                   1290:                        if ( length(MF) == 1 && deg(MF[0][0],TMP_S)==1 ) {
                   1291:                                P0 = cons(subst(MF[0][0],TMP_S,V[I]),P0);
                   1292:                        } else break;
                   1293:                }
                   1294:                if ( I == N ) {
                   1295:                        Indata = L[0]; AllData = L[1]; VData = L[2];
                   1296:                        GIN = Indata[0]; VDV = Indata[1]; WVDV = AllData[4];
                   1297:                        W = Indata[4];
                   1298:                        dp_set_weight(W);
                   1299:                        B = weyl_minipoly(GIN,VDV,0,WVDV);
                   1300:                        B = subst(B,s,-s-1);
                   1301:                        dp_set_weight(0);
                   1302:                        return [
                   1303:                          [P0,[1],cdr(fctr(B))],
                   1304:                          [[F],P0,[[s+1,1]]],
                   1305:                          [[0],[F],[]]
                   1306:                    ];
                   1307:                }
                   1308:        }
                   1309:
                   1310:        L2 = bf_strat_stage2(L);
                   1311:        S = bf_strat_stage3(L2);
                   1312:        R = [];
                   1313:        for ( T = S; T != []; T = cdr(T) ) {
                   1314:                Str = car(T);
                   1315:                B = Str[2];
                   1316:                B1 = [];
                   1317:                for ( U = B; U != []; U = cdr(U) )
                   1318:                        B1 = cons([subst(car(U)[0],s,-s-1),car(U)[1]],B1);
                   1319:                B1 = reverse(B1);
                   1320:                R = cons([Str[0],Str[1],B1],R);
                   1321:        }
                   1322:        return reverse(R);
                   1323: }
                   1324:
                   1325: /* returns [QQ,V,V0,B,BF,W] */
                   1326: /* QQ : ideal in C[x,s] (s=tdt), V=[x1,..,xn,t], V0 = [x1,..,xn] */
                   1327: /* B : global b-function, BF : factor list of B, W : weight */
                   1328:
                   1329: def bf_strat_stage2(L)
                   1330: {
                   1331:        T0 = time();
                   1332:        InData = L[0]; VData = L[2];
                   1333:        G1 = InData[0]; VDV = InData[1]; W = InData[4]; W0 = VData[4];
                   1334:        N = length(VDV); N1 = N/2;
                   1335:        V = InData[2]; DV = InData[3];
                   1336:        T = VData[2]; DT = VData[3];
                   1337:        V0 = VData[0]; DVR = VData[1];
                   1338:        dp_set_weight(W);
                   1339:        for ( I = 0; DVR != []; I++, DVR = cdr(DVR) ) {
                   1340:                DVRV = cons(DT,append(cdr(DVR),V));
                   1341:                M = elim_mat(VDV,DVRV);
                   1342:                for ( K = 0; K < N; K++ )
                   1343:                        M[1][K] = W[K];
1.2       noro     1344:                dp_ord(0); D1 = map(dp_ptod,G1,VDV);
                   1345:                H1 = map(dp_ht,D1); HC1 = map(dp_hc,D1);
1.1       noro     1346:                dp_ord(M); H2 = map(dp_ht,map(dp_ptod,G1,VDV));
                   1347:                if ( H1 == H2 )
                   1348:                        G2 = G1;
                   1349:                else
1.2       noro     1350:                        G2 = nd_gb_candidate(G1,VDV,M,0,HC1,1);
1.1       noro     1351:                G1 = elimination(G2,DVRV);
                   1352:        }
                   1353:        T1 = time();
                   1354:        B = weyl_minipoly(G1,VDV,0,T*DT);
                   1355:        T2 = time();
                   1356:        BF = cdr(fctr(B));
                   1357:
                   1358:        dp_set_weight(0);
                   1359:        G1 = map(psi0,G1,T,DT);
                   1360:        QQ = map(subst,map(b_subst,map(subst,G1,DT,1),T),T,var(B));
                   1361:        if ( type(getopt(ideal)) != -1 ) return [QQ,V];
                   1362:        print(["elim",(T1[0]+T1[1])-(T0[0]+T0[1])]);
                   1363:        print(["globalb",(T2[0]+T2[1])-(T1[0]+T1[1])]);
                   1364:        return [QQ,V,V0,B,BF,W0,DV];
                   1365: }
                   1366:
                   1367: def bf_strat_stage3(L)
                   1368: {
                   1369:        T0 = time();
                   1370:        QQ = L[0]; V0 = L[2]; B = L[3]; BF = L[4]; W0 = L[5];
                   1371:        NF = length(BF);
                   1372:        Data = vector(NF);
1.2       noro     1373:        W1 = W0? cons(1,append(W0,[1])) : 0;
1.1       noro     1374:        for ( I = J = 0; I < NF; I++ ) {
                   1375:                DI = tower_in_p(QQ,B,BF[I],V0,W0);
                   1376:                NDI = length(DI);
1.2       noro     1377:                dp_set_weight(W1);
1.1       noro     1378:                for ( K = 0; K < J; K++ ) {
                   1379:                        if ( length(DK=Data[K]) == NDI ) {
                   1380:                                for ( L = 0; L < NDI; L++ ) {
                   1381:                                        CL = DI[L][1]; CK = DK[L][1];
                   1382:                                        if ( !zero_inclusion(CL,CK,V0)
                   1383:                                                || !zero_inclusion(CK,CL,V0) ) break;
                   1384:                                }
                   1385:                                if ( L == NDI ) break;
                   1386:                        }
                   1387:                }
1.2       noro     1388:                dp_set_weight(0);
1.1       noro     1389:                if ( K < J ) {
                   1390:                        for ( L = 0, T = []; L < NDI; L++ )
                   1391:                                T = cons([[DK[L][0][0]*DI[L][0][0],DK[L][0][1]],
                   1392:                                        DK[L][1],DK[L][2]],T);
                   1393:                        Data[K] = reverse(T);
                   1394:                } else
                   1395:                        Data[J++] = DI;
                   1396:        }
                   1397:        Data1 = vector(J);
                   1398:        for ( I = 0; I < J; I++ )
                   1399:                Data1[I] = Data[I];
                   1400:        T1 = time();
1.2       noro     1401:        Str = stratify_bf(Data1,V0,W0);
1.1       noro     1402:        T2 = time();
                   1403:        print(["tower",(T1[0]+T1[1])-(T0[0]+T0[1])]);
                   1404:        print(["strat",(T2[0]+T2[1])-(T1[0]+T1[1])]);
                   1405:        return Str;
                   1406: }
                   1407:
                   1408: /*
                   1409:  InData = [GIN,VDV,V,DV,WtV]
                   1410:  AllData = [G0,GIN0,VDV0,W,WVDV,WtV0]
                   1411: */
                   1412:
                   1413: def bf_local(F,P)
                   1414: {
                   1415:        if ( type(Heu=getopt(heuristic)) == -1 ) Heu = 0;
                   1416:        if ( type(Vord=getopt(vord)) == -1 || type(Vord) != 4 ) Vord = 0;
                   1417:        if ( type(Wt=getopt(weight)) == -1 ) Wt = 0;
1.3       noro     1418:        if ( type(Op=getopt(op)) == -1 ) Op = 0;
1.1       noro     1419:        L = in_ww(F|weight=Wt,heuristic=Heu,vord=Vord);
                   1420:        InData = L[0]; AllData = L[1]; VData = L[2];
                   1421:        G = InData[0]; VDV = InData[1];
                   1422:        V = InData[2]; DV = InData[3];
                   1423:
                   1424:        V0 = VData[0]; DV0 = VData[1]; T = VData[2]; DT = VData[3]; W0 = VData[4];
                   1425:
                   1426:        L2 = bf_strat_stage2(L);
                   1427:
                   1428: /* L2 = [QQ,V,V0,B,BF,W] */
                   1429: /* QQ : ideal in C[x,s] (s=tdt), V=[t,x1,..,xn], V0 = [x1,..,xn] */
                   1430: /* B : global b-function, BF : factor list of B, W : weight */
                   1431:
                   1432:        QQ = L2[0]; B = L2[3]; BF = L2[4]; W = L2[5];
                   1433:
                   1434:        NF = length(BF);
                   1435:        BP = [];
                   1436:        dp_set_weight(0);
                   1437:        for ( I = J = 0; I < NF; I++ ) {
                   1438:                List = compute_exponent(QQ,B,BF[I],P,V0,W0);
                   1439:                DI = List[0]; QQI = List[1];
                   1440:                if ( DI )
                   1441:                        BP = cons([BF[I][0],DI],BP);
                   1442:                if ( I == 0 )
                   1443:                        Id = QQI;
                   1444:                else
                   1445:                        Id = ideal_intersection(Id,QQI,V0,0);
                   1446:        }
                   1447:        for ( List = Id; List != []; List = cdr(List) )
                   1448:                if ( subst_vars(car(List),P) )
                   1449:                        break;
                   1450:        if ( List == [] ) error("bf_local : inconsitent intersection");
                   1451:        Ax = car(List);
1.3       noro     1452:        LB = [];
                   1453:        for ( BPT = 1, List = BP; List != []; List = cdr(List) ) {
1.1       noro     1454:                BPT *= car(List)[0]^car(List)[1];
1.3       noro     1455:                LB = cons([subst(car(List)[0],s,-s-1),car(List)[1]],LB);
                   1456:        }
                   1457:        LB = reverse(LB);
                   1458:        if ( !Op ) return LB;
                   1459:
1.1       noro     1460:        BPT = weyl_subst(BPT,T*DT,VDV);
                   1461:
                   1462:        /* computation using G0,GIN0,VDV0 */
                   1463:        G0 = AllData[0]; GIN0 = AllData[1]; VDV0 = AllData[2]; WtV0 = AllData[5];
                   1464:        dp_set_weight(WtV0); dp_ord(0);
                   1465:        PS = map(dp_ptod,GIN0,VDV0); Len = length(PS);
                   1466:        for ( I = Len-1, Ind = []; I >= 0; I-- ) Ind = cons(I,Ind);
                   1467:        /* QR = [D,M,Coef] */
                   1468:        AxBPT = dp_ptod(Ax*BPT,VDV0);
                   1469:        QR = weyl_nf_quo(Ind,AxBPT,1,PS);
                   1470:        if ( !weyl_nf_quo_check(AxBPT,PS,QR) ) error("bf_local : invalid quotient");
                   1471:        if ( QR[0] ) error("bf_local : invalid quotient");
                   1472:        Den = QR[1]; Coef = QR[2];
                   1473:        for ( I = 0, R = Den*AxBPT; I < Len; I++ )
                   1474:                R -= dp_weyl_mul(Coef[I],dp_ptod(G0[I],VDV0));
                   1475:        R = dp_dtop(R,VDV0);
                   1476:        CR = conv_tdt(R,F,V0,DV0,T,DT);
                   1477:
                   1478:        dp_set_weight(0);
1.3       noro     1479:        return [LB,Ax,CR];
1.1       noro     1480: }
                   1481:
                   1482: /* t^(l+k)*dt^l (k>l) -> (s-k)(s-k-1)...(s-(k+l-1))t^k */
                   1483: def conv_tdt(P,F,V0,DV0,T,DT)
                   1484: {
                   1485:        DP = dp_ptod(P,[T,DT]);
                   1486:        VDV = append(cons(T,V0),cons(DT,DV0));
                   1487:        R = 0;
                   1488:        DF = dp_ptod(F,VDV);
                   1489:        for ( ; DP; DP = dp_rest(DP) ) {
                   1490:                C = dp_hc(DP);
                   1491:                E = dp_etov(dp_ht(DP));
                   1492:                L = E[1]; K = E[0]-E[1];
                   1493:                S = 1;
                   1494:                for ( I = 0; I < L; I++ )
                   1495:                        S *= ((-T-1)-K-I);
                   1496:                U = dp_ptod(C*S,VDV);
                   1497:                for ( I = 1; I < K; I++ )
                   1498:                        U = dp_weyl_mul(U,DF);
                   1499:                R += dp_dtop(U,VDV);
                   1500:        }
                   1501:        return subst(R,T,s);
                   1502: }
                   1503:
1.2       noro     1504: /* W1=[W,1], W2=[1,W,1] */
                   1505:
                   1506: def merge_tower(Str,Tower,V,W1,W2)
1.1       noro     1507: {
                   1508:        Prev = car(Tower); T = cdr(Tower);
                   1509:        Str1 = [];
                   1510:        for ( ; T != []; T = cdr(T) ) {
                   1511:                Cur = car(T);
                   1512:                Str1 = cons([Cur[1],Prev[1],[Prev[0]]],Str1);
                   1513:                Prev = Cur;
                   1514:        }
                   1515:        Str1 = cons([[0],Prev[1],[]],Str1);
                   1516:        Str1 = reverse(Str1);
                   1517:        if ( Str == [] ) return Str1;
                   1518:        U = [];
                   1519:        for ( T = Str; T != []; T = cdr(T) ) {
                   1520:                for ( S = Str1; S != []; S = cdr(S) ) {
1.2       noro     1521:                        Int = int_cons(T[0],S[0],V,W1,W2);
1.1       noro     1522:                        if ( Int[0] != [1] )
                   1523:                                U = cons(append(Int,[append(T[0][2],S[0][2])]),U);
                   1524:                }
                   1525:        }
                   1526:        U = reverse(U);
                   1527:        return U;
                   1528: }
                   1529:
1.2       noro     1530: def stratify_bf(Data,V,W)
1.1       noro     1531: {
                   1532:        N = length(Data);
                   1533:        R = [];
1.2       noro     1534:        if ( W ) {
                   1535:                W1 = append(W,[1]);
                   1536:                W2 = cons(1,W1);
                   1537:        } else
                   1538:                W1 = W2 = 0;
1.1       noro     1539:        for ( I = 0; I < N; I++ )
1.2       noro     1540:                R = merge_tower(R,Data[I],V,W1,W2);
1.1       noro     1541:        return R;
                   1542: }
                   1543:
                   1544: def tdt_homogenize(F,L)
                   1545: {
                   1546:        TY1 = L[0]; T = TY1[0]; Y1 = TY1[1];
                   1547:        TY2 = L[1]; DT = TY2[0]; Y2 = TY2[1];
                   1548:        DF = dp_ptod(F,[T,DT]);
                   1549:        for ( R = 0; DF; DF = dp_rest(DF) ) {
                   1550:                M = dp_hm(DF);
                   1551:                E = dp_etov(M);
                   1552:                W = E[1]-E[0];
                   1553:                if ( W > 0 ) R += Y1^W*dp_dtop(M,[T,DT]);
                   1554:                else R += Y2^W*dp_dtop(M,[T,DT]);
                   1555:        }
                   1556:        return R;
                   1557: }
                   1558:
                   1559: def sing(F,V)
                   1560: {
                   1561:        R = [F];
                   1562:        for ( T = V; T != []; T = cdr(T) )
                   1563:                R = cons(diff(F,car(T)),R);
                   1564:        G = nd_gr_trace(R,V,1,1,0);
                   1565:        return G;
                   1566: }
                   1567:
                   1568: def tower_in_p(B,F,FD,V,W)
                   1569: {
                   1570:        TT = ttttt;
                   1571:        N = length(V); S = var(F); SV = cons(S,V); V1 = cons(TT,SV);
                   1572:        Wt = append(append([1,1],W),[1]);
                   1573:        dp_set_weight(Wt);
                   1574:
                   1575:        F1 = FD[0]; D = FD[1];
                   1576:        O1 = [[0,1],[0,N+1]]; O2 = [[0,1],[0,N]];
                   1577:
                   1578:        TF = sdiv(F,F1^D);
                   1579:
                   1580:        T = nd_gr_trace(cons((1-TT)*TF,vtol(TT*ltov(B))),V1,1,1,O1);
                   1581:        T = elimination(T,SV);
                   1582:        Q = map(sdiv,T,TF);
                   1583:        dp_set_weight(cdr(Wt));
                   1584:        QQ = elimination(nd_gr(Q,SV,0,O2),V);
                   1585:        E = [[[F1,0],QQ,PD]];
                   1586:
                   1587:        for ( I = D-1; I >= 0; I-- ) {
                   1588:            dp_set_weight(Wt);
                   1589:                T = nd_gr_trace(cons((1-TT)*F1,vtol(TT*ltov(Q))),V1,1,1,O1);
                   1590:                T = elimination(T,SV);
                   1591:                Q = map(sdiv,T,F1);
                   1592:            dp_set_weight(cdr(Wt));
                   1593:                QQ = elimination(nd_gr(Q,SV,0,O2),V);
                   1594:                E = cons([[F1,D-I],QQ,PD],E);
                   1595:        }
                   1596:        dp_set_weight(0);
                   1597:        return E;
                   1598: }
                   1599:
                   1600: def subst_vars(F,P)
                   1601: {
                   1602:        for ( ; P != []; P = cdr(cdr(P)) )
                   1603:                F = subst(F,P[0],P[1]);
                   1604:        return F;
                   1605: }
                   1606:
                   1607: def compute_exponent(B,F,FD,P,V,W)
                   1608: {
                   1609:        TT = ttttt;
                   1610:        N = length(V); S = var(F); SV = cons(S,V); V1 = cons(TT,SV);
                   1611:        F1 = FD[0]; D = FD[1];
                   1612:
                   1613:        Wt = append(append([1,1],W),[1]);
                   1614:        dp_set_weight(Wt);
                   1615:        O1 = [[0,1],[0,N+1]]; O2 = [[0,1],[0,N]];
                   1616:
                   1617:        TF = sdiv(F,F1^D);
                   1618:
                   1619:        T = nd_gr_trace(cons((1-TT)*TF,vtol(TT*ltov(B))),V1,0,1,O1);
                   1620:        T = elimination(T,SV);
                   1621:        Q = map(sdiv,T,TF);
                   1622:        QQ = elimination(nd_gr(Q,SV,0,O2),V);
                   1623:
                   1624:        for ( I = 0; I < D; I++ ) {
                   1625:                for ( T = QQ; T != []; T = cdr(T) )
                   1626:                        if ( subst_vars(car(T),P) ) {
                   1627:                                dp_set_weight(0);
                   1628:                                return [I,QQ];
                   1629:                        }
                   1630:                T = nd_gr_trace(cons((1-TT)*F1,vtol(TT*ltov(Q))),V1,0,1,O1);
                   1631:                T = elimination(T,SV);
                   1632:                Q = map(sdiv,T,F1);
                   1633:                QQ = elimination(nd_gr(Q,SV,0,O2),V);
                   1634:        }
                   1635:        dp_set_weight(0);
                   1636:        return [D,QQ];
                   1637: }
                   1638:
                   1639: /* V(B) subset V(A) ? */
                   1640:
                   1641: def zero_inclusion(A,B,V)
                   1642: {
                   1643:        NV = ttttt;
                   1644:        for ( T = A; T != []; T = cdr(T) ) {
                   1645:                F = car(T);
1.2       noro     1646:                G = nd_gr_trace(cons(NV*F-1,B),cons(NV,V),1,1,0);
1.1       noro     1647:                if ( type(car(G)) != 1 ) return 0;
                   1648:        }
                   1649:        return 1;
                   1650: }
                   1651:
                   1652: def weyl_divide_by_right(G,H,V,O)
                   1653: {
                   1654:        dp_ord(O); G = dp_ptod(G,V); H = dp_ptod(H,V);
                   1655:        CH = dp_hc(H);
                   1656:        for ( Q = 0; G; ) {
                   1657:                if ( !dp_redble(G,H) ) return 0;
                   1658:                CG = dp_hc(G);
                   1659:                Q1 = CG/CH*dp_subd(G,H);
                   1660:                G -= dp_weyl_mul(Q1,H);
                   1661:                Q += Q1;
                   1662:        }
                   1663:        return dp_dtop(Q,V);
                   1664: }
                   1665:
                   1666: def elim_mat(V,W)
                   1667: {
                   1668:        N = length(V);
                   1669:        M = matrix(N+1,N);
                   1670:        for ( J = 0; J < N; J++ ) if ( !member(V[J],W) ) M[0][J] = 1;
                   1671:        for ( J = 0; J < N; J++ ) M[1][J] = 1;
                   1672:        for ( I = 2; I <= N; I++ ) M[I][N-I+1] = -1;
                   1673:        return M;
                   1674: }
                   1675:
                   1676: /* (P-Q)cap(R-S)=(P cap Q^c)cap(R cap S^c)=(P cap R)cap(Q cup S)^c
                   1677:    =(P cap R)-(Q cup S)
                   1678: */
                   1679:
1.2       noro     1680: def int_cons(A,B,V,W1,W2)
1.1       noro     1681: {
                   1682:        AZ = A[0]; AN = A[1];
                   1683:        BZ = B[0]; BN = B[1];
1.2       noro     1684:        if ( W1 ) dp_set_weight(W1);
1.1       noro     1685:        CZ = nd_gr_trace(append(AZ,BZ),V,1,1,0);
1.2       noro     1686:        if ( W2 ) dp_set_weight(W2);
1.1       noro     1687:        CN = ideal_intersection(AN,BN,V,0);
1.2       noro     1688:        ZI = zero_inclusion(CN,CZ,V);
                   1689:        dp_set_weight(0);
                   1690:        if ( ZI )
1.1       noro     1691:                return [[1],[]];
                   1692:        else
                   1693:                return [CZ,CN];
                   1694: }
                   1695:
                   1696: def ideal_intersection(A,B,V,Ord)
                   1697: {
                   1698:        T = ttttt;
                   1699:        G = nd_gr_trace(append(vtol(ltov(A)*T),vtol(ltov(B)*(1-T))),
                   1700:                cons(T,V),1,1,[[0,1],[Ord,length(V)]]);
                   1701:        return elimination(G,V);
                   1702: }
1.2       noro     1703:
                   1704: def nd_gb_candidate(G,V,Ord,Homo,HC,Weyl)
                   1705: {
                   1706:        Ind = 0;
                   1707:        N = length(HC);
                   1708:        while ( 1 ) {
                   1709:                P = lprime(Ind++);
                   1710:                for ( I = 0; I < N; I++ )
                   1711:                        if ( !(HC[I]%P) ) break;
                   1712:                if ( I < N ) continue;
                   1713:                if ( Weyl )
                   1714:                        G = nd_weyl_gr_trace(G,V,Homo,-P,Ord);
                   1715:                else
                   1716:                        G = nd_gr_trace(G,V,Homo,-P,Ord);
                   1717:                if ( G ) return G;
                   1718:        }
                   1719: }
                   1720:
1.1       noro     1721: endmodule $
                   1722: end$
                   1723:
                   1724:

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