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

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

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