=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/gr.c,v retrieving revision 1.39 retrieving revision 1.65 diff -u -p -r1.39 -r1.65 --- OpenXM_contrib2/asir2000/builtin/gr.c 2001/11/19 00:57:10 1.39 +++ OpenXM_contrib2/asir2000/builtin/gr.c 2009/02/09 10:21:29 1.65 @@ -45,7 +45,7 @@ * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. * - * $OpenXM: OpenXM_contrib2/asir2000/builtin/gr.c,v 1.38 2001/11/16 10:35:07 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/builtin/gr.c,v 1.64 2007/09/19 05:56:01 noro Exp $ */ #include "ca.h" #include "parse.h" @@ -69,7 +69,7 @@ double get_rtime(); struct oEGT eg_nf,eg_nfm; struct oEGT eg_znfm,eg_pz,eg_np,eg_ra,eg_mc,eg_gc; -int TP,NBP,NMP,NFP,NDP,ZR,NZR; +int TP,N_BP,NMP,NFP,NDP,ZR,NZR; extern int (*cmpdl)(); extern int do_weyl; @@ -89,26 +89,28 @@ static VL VC; int PCoeffs; int DP_Print = 0; +int DP_PrintShort = 0; int DP_Multiple = 0; int DP_NFStat = 0; LIST Dist = 0; int NoGCD = 0; int GenTrace = 0; +int GenSyz = 0; int OXCheck = -1; -static int NoSugar = 0; +int NoSugar = 0; static int NoCriB = 0; static int NoGC = 0; static int NoMC = 0; static int NoRA = 0; -static int DP_PrintShort = 0; static int ShowMag = 0; static int Stat = 0; -static int Denominator = 1; -static int Top = 0; -static int Reverse = 0; +int Denominator = 1; +int Top = 0; +int Reverse = 0; static int Max_mag = 0; -static char *Demand = 0; +static int Max_coef = 0; +char *Demand = 0; static int PtozpRA = 0; int doing_f4; @@ -122,17 +124,70 @@ void Pox_pop_local(NODE,Obj *); INLINE int eqdl(int nv,DL dl1,DL dl2) { int i; - int *b1,*b2; + int *p1,*p2; if ( dl1->td != dl2->td ) return 0; - for ( i = 0, b1 = dl1->d, b2 = dl2->d; i < nv; i++ ) - if ( b1[i] != b2[i] ) - break; - if ( i == nv ) - return 1; - else - return 0; + i = nv-1; + p1 = dl1->d; + p2 = dl2->d; + while ( i >= 7 ) { + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + i -= 8; + } + switch ( i ) { + case 6: + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + return 1; + case 5: + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + return 1; + case 4: + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + return 1; + case 3: + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + return 1; + case 2: + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + return 1; + case 1: + if ( *p1++ != *p2++ ) return 0; + if ( *p1++ != *p2++ ) return 0; + return 1; + case 0: + if ( *p1++ != *p2++ ) return 0; + return 1; + default: + return 1; + } } /* b[] should be cleared */ @@ -244,12 +299,13 @@ void pdl(NODE f) void dp_gr_main(LIST f,LIST v,Num homo,int modular,int field,struct order_spec *ord,LIST *rp) { int i,mindex,m,nochk; - struct order_spec ord1; + struct order_spec *ord1; Q q; VL fv,vv,vc; NODE fd,fd0,fi,fi0,r,r0,t,subst,x,s,xx; NODE ind,ind0; LIST trace,gbindex; + int input_is_dp = 0; mindex = 0; nochk = 0; dp_fcoeffs = field; get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc); @@ -262,14 +318,22 @@ void dp_gr_main(LIST f,LIST v,Num homo,int modular,int homogenize_order(ord,NVars,&ord1); for ( fd0 = fi0 = 0, t = BDY(f); t; t = NEXT(t) ) { NEXTNODE(fd0,fd); NEXTNODE(fi0,fi); - ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fi)); dp_homo((DP)BDY(fi),(DP *)&BDY(fd)); + if ( BDY(t) && OID(BDY(t)) == O_DP ) { + dp_sort((DP)BDY(t),(DP *)&BDY(fi)); input_is_dp = 1; + } else + ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fi)); + dp_homo((DP)BDY(fi),(DP *)&BDY(fd)); } if ( fd0 ) NEXT(fd) = 0; if ( fi0 ) NEXT(fi) = 0; - initd(&ord1); + initd(ord1); } else { for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) { - NEXTNODE(fd0,fd); ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd)); + NEXTNODE(fd0,fd); + if ( BDY(t) && OID(BDY(t)) == O_DP ) { + dp_sort((DP)BDY(t),(DP *)&BDY(fd)); input_is_dp = 1; + } else + ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd)); } if ( fd0 ) NEXT(fd) = 0; fi0 = fd0; @@ -286,7 +350,7 @@ void dp_gr_main(LIST f,LIST v,Num homo,int modular,int init_stat(); while ( 1 ) { if ( homo ) { - initd(&ord1); CNVars = NVars+1; + initd(ord1); CNVars = NVars+1; } if ( DP_Print && modular ) { fprintf(asir_out,"mod= %d, eval = ",m); printsubst(subst); @@ -317,7 +381,10 @@ void dp_gr_main(LIST f,LIST v,Num homo,int modular,int } for ( r0 = 0, ind0 = 0; x; x = NEXT(x) ) { NEXTNODE(r0,r); dp_load((int)BDY(x),&ps[(int)BDY(x)]); - dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r)); + if ( input_is_dp ) + BDY(r) = (pointer)ps[(int)BDY(x)]; + else + dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r)); NEXTNODE(ind0,ind); STOQ((int)BDY(x),q); BDY(ind) = q; } @@ -338,16 +405,63 @@ void dp_gr_main(LIST f,LIST v,Num homo,int modular,int } print_stat(); if ( ShowMag ) - fprintf(asir_out,"\nMax_mag=%d\n",Max_mag); + fprintf(asir_out,"\nMax_mag=%d, Max_coef=%d\n",Max_mag, Max_coef); } +void dp_interreduce(LIST f,LIST v,int field,struct order_spec *ord,LIST *rp) +{ + int i,mindex,m,nochk; + struct order_spec *ord1; + Q q; + VL fv,vv,vc; + NODE fd,fd0,fi,fi0,r,r0,t,subst,x,s,xx; + NODE ind,ind0; + LIST trace,gbindex; + int input_is_dp = 0; + + mindex = 0; nochk = 0; dp_fcoeffs = field; + get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc); + NVars = length((NODE)vv); PCoeffs = vc ? 1 : 0; VC = vc; + CNVars = NVars; + if ( ord->id && NVars != ord->nv ) + error("dp_interreduce : invalid order specification"); + initd(ord); + for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) { + NEXTNODE(fd0,fd); + if ( BDY(t) && OID(BDY(t)) == O_DP ) { + dp_sort((DP)BDY(t),(DP *)&BDY(fd)); input_is_dp = 1; + } else + ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd)); + } + if ( fd0 ) NEXT(fd) = 0; + fi0 = fd0; + setup_arrays(fd0,0,&s); + init_stat(); + x = s; + reduceall(x,&xx); x = xx; + for ( r0 = 0, ind0 = 0; x; x = NEXT(x) ) { + NEXTNODE(r0,r); dp_load((int)BDY(x),&ps[(int)BDY(x)]); + if ( input_is_dp ) + BDY(r) = (pointer)ps[(int)BDY(x)]; + else + dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r)); + NEXTNODE(ind0,ind); + STOQ((int)BDY(x),q); BDY(ind) = q; + } + if ( r0 ) NEXT(r) = 0; + if ( ind0 ) NEXT(ind) = 0; + MKLIST(*rp,r0); + MKLIST(gbindex,ind0); +} + void dp_gr_mod_main(LIST f,LIST v,Num homo,int m,struct order_spec *ord,LIST *rp) { - struct order_spec ord1; + struct order_spec *ord1; VL fv,vv,vc; NODE fd,fd0,r,r0,t,x,s,xx; DP a,b,c; -extern struct oEGT eg_red_mod; + extern struct oEGT eg_red_mod; + int input_is_dp = 0; get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc); NVars = length((NODE)vv); PCoeffs = vc ? 1 : 0; VC = vc; @@ -357,7 +471,11 @@ extern struct oEGT eg_red_mod; initd(ord); if ( homo ) { for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) { - ptod(CO,vv,(P)BDY(t),&a); dp_homo(a,&b); + if ( BDY(t) && OID(BDY(t)) == O_DP ) { + dp_sort((DP)BDY(t),&a); input_is_dp = 1; + } else + ptod(CO,vv,(P)BDY(t),&a); + dp_homo(a,&b); if ( PCoeffs ) dp_mod(b,m,0,&c); else @@ -366,10 +484,13 @@ extern struct oEGT eg_red_mod; NEXTNODE(fd0,fd); BDY(fd) = (pointer)c; } } - homogenize_order(ord,NVars,&ord1); initd(&ord1); + homogenize_order(ord,NVars,&ord1); initd(ord1); } else { for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) { - ptod(CO,vv,(P)BDY(t),&b); + if ( BDY(t) && OID(BDY(t)) == O_DP ) { + dp_sort((DP)BDY(t),&b); input_is_dp = 1; + } else + ptod(CO,vv,(P)BDY(t),&b); if ( PCoeffs ) dp_mod(b,m,0,&c); else @@ -383,7 +504,7 @@ extern struct oEGT eg_red_mod; setup_arrays(fd0,m,&s); init_stat(); if ( homo ) { - initd(&ord1); CNVars = NVars+1; + initd(ord1); CNVars = NVars+1; } /* init_eg(&eg_red_mod); */ x = gb_mod(s,m); @@ -395,11 +516,19 @@ extern struct oEGT eg_red_mod; reduceall_mod(x,m,&xx); x = xx; if ( PCoeffs ) for ( r0 = 0; x; x = NEXT(x) ) { - NEXTNODE(r0,r); mdtop(CO,m,vv,ps[(int)BDY(x)],(P *)&BDY(r)); + NEXTNODE(r0,r); + if ( input_is_dp ) + mdtodp(ps[(int)BDY(x)],(DP *)&BDY(r)); + else + mdtop(CO,m,vv,ps[(int)BDY(x)],(P *)&BDY(r)); } else for ( r0 = 0; x; x = NEXT(x) ) { - NEXTNODE(r0,r); _dtop_mod(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r)); + NEXTNODE(r0,r); + if ( input_is_dp ) + _mdtodp(ps[(int)BDY(x)],(DP *)&BDY(r)); + else + _dtop_mod(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r)); } print_stat(); if ( r0 ) NEXT(r) = 0; @@ -411,6 +540,7 @@ void dp_f4_main(LIST f,LIST v,struct order_spec *ord,L int homogen; VL fv,vv,vc; NODE fd,fd0,r,r0,t,x,s,xx; + int input_is_dp = 0; dp_fcoeffs = 0; get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc); @@ -420,7 +550,11 @@ void dp_f4_main(LIST f,LIST v,struct order_spec *ord,L error("dp_f4_main : invalid order specification"); initd(ord); for ( fd0 = 0, t = BDY(f), homogen = 1; t; t = NEXT(t) ) { - NEXTNODE(fd0,fd); ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd)); + NEXTNODE(fd0,fd); + if ( BDY(t) && OID(BDY(t)) == O_DP ) { + dp_sort((DP)BDY(t),(DP *)&BDY(fd)); input_is_dp = 1; + } else + ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd)); if ( homogen ) homogen = dp_homogeneous(BDY(fd)); } @@ -432,7 +566,10 @@ void dp_f4_main(LIST f,LIST v,struct order_spec *ord,L } for ( r0 = 0; x; x = NEXT(x) ) { NEXTNODE(r0,r); dp_load((int)BDY(x),&ps[(int)BDY(x)]); - dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r)); + if ( input_is_dp ) + BDY(r) = (pointer)ps[(int)BDY(x)]; + else + dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r)); } if ( r0 ) NEXT(r) = 0; MKLIST(*rp,r0); @@ -444,6 +581,7 @@ void dp_f4_mod_main(LIST f,LIST v,int m,struct order_s VL fv,vv,vc; DP b,c,c1; NODE fd,fd0,r,r0,t,x,s,xx; + int input_is_dp = 0; dp_fcoeffs = 0; get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc); @@ -453,7 +591,10 @@ void dp_f4_mod_main(LIST f,LIST v,int m,struct order_s error("dp_f4_mod_main : invalid order specification"); initd(ord); for ( fd0 = 0, t = BDY(f), homogen = 1; t; t = NEXT(t) ) { - ptod(CO,vv,(P)BDY(t),&b); + if ( BDY(t) && OID(BDY(t)) == O_DP ) { + dp_sort((DP)BDY(t),&b); input_is_dp = 1; + } else + ptod(CO,vv,(P)BDY(t),&b); if ( homogen ) homogen = dp_homogeneous(b); _dp_mod(b,m,0,&c); @@ -473,7 +614,11 @@ void dp_f4_mod_main(LIST f,LIST v,int m,struct order_s reduceall_mod(x,m,&xx); x = xx; } for ( r0 = 0; x; x = NEXT(x) ) { - NEXTNODE(r0,r); _dtop_mod(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r)); + NEXTNODE(r0,r); + if ( input_is_dp ) + _mdtodp(ps[(int)BDY(x)],(DP *)&BDY(r)); + else + _dtop_mod(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r)); } if ( r0 ) NEXT(r) = 0; MKLIST(*rp,r0); @@ -600,6 +745,8 @@ NODE gb_f4(NODE f) /* initial bases are monic */ unsigned int **psca; +GeoBucket create_bucket(); +DL remove_head_bucket(GeoBucket,int); NODE gb_f4_mod(NODE f,int m) { @@ -618,14 +765,16 @@ NODE gb_f4_mod(NODE f,int m) int *indred,*isred; CDP ri; int pscalen; - struct oEGT tmp0,tmp1,eg_split_symb,eg_split_elim1,eg_split_elim2; - extern struct oEGT eg_symb,eg_elim1,eg_elim2; + GeoBucket bucket; + DL head; + struct oEGT tmp0,tmp1,eg_split_symb,eg_split_conv,eg_split_elim1,eg_split_elim2; + extern struct oEGT eg_symb,eg_conv,eg_elim1,eg_elim2; /* initialize coeffcient array list of ps[] */ pscalen = pslen; psca = (unsigned int **)MALLOC(pscalen*sizeof(unsigned int *)); - init_eg(&eg_symb); init_eg(&eg_elim1); init_eg(&eg_elim2); + init_eg(&eg_symb); init_eg(&eg_conv); init_eg(&eg_elim1); init_eg(&eg_elim2); for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) { i = (int)BDY(r); d = updpairs(d,g,i); @@ -640,41 +789,57 @@ NODE gb_f4_mod(NODE f,int m) minsugar(d,&dm,&dr); d = dr; if ( DP_Print ) fprintf(asir_out,"sugar=%d\n",dm->sugar); - blist = 0; s0 = 0; + blist = 0; + bucket = create_bucket(); /* asph : sum of all head terms of spoly */ for ( t = dm; t; t = NEXT(t) ) { _dp_sp_mod(ps[t->dp1],ps[t->dp2],m,&sp); /* fprintf(stderr,"splen=%d-",dp_nt(sp)); */ if ( sp ) { MKNODE(bt,sp,blist); blist = bt; - s0 = symb_merge(s0,dp_dllist(sp),nv); + add_bucket(bucket,dp_dllist(sp),nv); /* fprintf(stderr,"%d-",length(s0)); */ } } +#if 0 if ( DP_Print ) fprintf(asir_out,"initial spmat : %d x %d ",length(blist),length(s0)); +#endif /* s0 : all the terms appeared in symbolic reduction */ - for ( s = s0, nred = 0; s; s = NEXT(s) ) { + nred = 0; + s0 = 0; + while ( 1 ) { + head = remove_head_bucket(bucket,nv); + if ( !head ) break; + else { + NEXTNODE(s0,s); + BDY(s) = (pointer)head; + } for ( r = gall; r; r = NEXT(r) ) - if ( _dl_redble(BDY(ps[(int)BDY(r)])->dl,BDY(s),nv) ) + if ( _dl_redble(BDY(ps[(int)BDY(r)])->dl,head,nv) ) break; if ( r ) { - dltod(BDY(s),nv,&tdp); + dltod(head,nv,&tdp); dp_subd(tdp,ps[(int)BDY(r)],&sd); dt = mul_dllist(BDY(sd)->dl,ps[(int)BDY(r)]); + add_bucket(bucket,NEXT(dt),nv); /* fprintf(stderr,"[%d]",length(dt)); */ /* list of [t,f] */ bt1 = mknode(2,BDY(sd)->dl,BDY(r)); MKNODE(bt,bt1,blist); blist = bt; - symb_merge(s,dt,nv); -/* fprintf(stderr,"%d-",length(s)); */ +/* fprintf(stderr,"%d-",length(s0)); */ nred++; } } + if ( s0 ) NEXT(s) = 0; /* fprintf(stderr,"\n"); */ + get_eg(&tmp1); add_eg(&eg_symb,&tmp0,&tmp1); + init_eg(&eg_split_symb); add_eg(&eg_split_symb,&tmp0,&tmp1); + if ( DP_Print ) fprintf(asir_out,"number of reducers : %d\n",nred); + get_eg(&tmp0); /* the first nred polys in blist are reducers */ /* row = the number of all the polys */ for ( r = blist, row = 0; r; r = NEXT(r), row++ ); @@ -713,8 +878,8 @@ NODE gb_f4_mod(NODE f,int m) for ( j = 0, k = 0; j < col; j++ ) if ( !isred[j] ) st[k++] = at[j]; - get_eg(&tmp1); add_eg(&eg_symb,&tmp0,&tmp1); - init_eg(&eg_split_symb); add_eg(&eg_split_symb,&tmp0,&tmp1); + get_eg(&tmp1); add_eg(&eg_conv,&tmp0,&tmp1); + init_eg(&eg_split_conv); add_eg(&eg_split_conv,&tmp0,&tmp1); get_eg(&tmp1); /* spoly matrix; stored in reduced form; terms in ht[] are omitted */ @@ -771,6 +936,7 @@ NODE gb_f4_mod(NODE f,int m) if ( DP_Print ) { fprintf(asir_out,"done rank = %d\n",rank,row,col); print_eg("Symb",&eg_split_symb); + print_eg("Conv",&eg_split_conv); print_eg("Elim1",&eg_split_elim1); print_eg("Elim2",&eg_split_elim2); fprintf(asir_out,"\n"); @@ -813,6 +979,7 @@ NODE gb_f4_mod(NODE f,int m) } if ( DP_Print ) { print_eg("Symb",&eg_symb); + print_eg("Conv",&eg_conv); print_eg("Elim1",&eg_elim1); print_eg("Elim2",&eg_elim2); fflush(asir_out); @@ -1019,6 +1186,20 @@ void pltovl(LIST l,VL *vl) *vl = r0; } +void vltopl(VL vl,LIST *l) +{ + VL n; + NODE r,r0; + P p; + + n = vl; + for ( r0 = 0; n; n = NEXT(n) ) { + NEXTNODE(r0,r); MKV(n->v,p); BDY(r) = (pointer)p; + } + if ( r0 ) NEXT(r) = 0; + MKLIST(*l,r0); +} + void makesubst(VL v,NODE *s) { NODE r,r0; @@ -1158,11 +1339,9 @@ void prim_part(DP f,int m,DP *r) if ( PCoeffs ) dp_prim_mod(f,m,r); else - *r = f; + _dp_monic(f,m,r); } else { - if ( dp_fcoeffs ) - *r = f; - else if ( PCoeffs ) + if ( dp_fcoeffs || PCoeffs ) dp_prim(f,r); else dp_ptozp(f,r); @@ -1251,9 +1430,10 @@ void reduceall(NODE in,NODE *h) w = (int *)ALLOCA(n*sizeof(int)); for ( i = 0, t = r; i < n; i++, t = NEXT(t) ) w[i] = (int)BDY(t); + /* w[i] < 0 : reduced to 0 */ for ( i = 0; i < n; i++ ) { for ( top = 0, j = n-1; j >= 0; j-- ) - if ( j != i ) { + if ( w[j] >= 0 && j != i ) { MKNODE(t,(pointer)w[j],top); top = t; } get_eg(&tmp0); @@ -1275,10 +1455,16 @@ void reduceall(NODE in,NODE *h) if ( DP_Print || DP_PrintShort ) { fprintf(asir_out,"."); fflush(asir_out); } - w[i] = newps(g1,0,(NODE)0); + if ( g1 ) { + w[i] = newps(g1,0,(NODE)0); + } else { + w[i] = -1; + } } for ( top = 0, j = n-1; j >= 0; j-- ) { - MKNODE(t,(pointer)w[j],top); top = t; + if ( w[j] >= 0 ) { + MKNODE(t,(pointer)w[j],top); top = t; + } } *h = top; if ( DP_Print || DP_PrintShort ) @@ -1304,9 +1490,10 @@ void reduceall_mod(NODE in,int m,NODE *h) w = (int *)ALLOCA(n*sizeof(int)); for ( i = 0, t = r; i < n; i++, t = NEXT(t) ) w[i] = (int)BDY(t); + /* w[i] < 0 : reduced to 0 */ for ( i = 0; i < n; i++ ) { for ( top = 0, j = n-1; j >= 0; j-- ) - if ( j != i ) { + if ( w[j] >= 0 && j != i ) { MKNODE(t,(pointer)w[j],top); top = t; } get_eg(&tmp0); @@ -1320,10 +1507,16 @@ void reduceall_mod(NODE in,int m,NODE *h) if ( DP_Print || DP_PrintShort ) { fprintf(asir_out,"."); fflush(asir_out); } - w[i] = newps_mod(g,m); + if ( g ) { + w[i] = newps_mod(g,m); + } else { + w[i] = -1; + } } for ( top = 0, j = n-1; j >= 0; j-- ) { - MKNODE(t,(pointer)w[j],top); top = t; + if ( w[j] >= 0 ) { + MKNODE(t,(pointer)w[j],top); top = t; + } } *h = top; if ( DP_Print || DP_PrintShort ) @@ -1537,7 +1730,7 @@ void minsugar(DP_pairs d,DP_pairs *dm,DP_pairs *dr) NODE gb(NODE f,int m,NODE subst) { - int i,nh,prev,mag; + int i,nh,prev,mag,mag0,magt; NODE r,g,gall; DP_pairs d; DP_pairs l; @@ -1551,8 +1744,10 @@ NODE gb(NODE f,int m,NODE subst) static prev_sugar = -1; Max_mag = 0; + Max_coef = 0; prev = 1; doing_f4 = 0; + init_denomlist(); if ( m ) { psm = (DP *)MALLOC(pslen*sizeof(DP)); for ( i = 0; i < psn; i++ ) @@ -1602,7 +1797,7 @@ NODE gb(NODE f,int m,NODE subst) _dp_nf(gall,h,ps,!Top,&nf); else _dp_nf_z(gall,h,ps,!Top,DP_Multiple,&nf); - if ( DP_Print ) + if ( DP_Print && nf ) fprintf(asir_out,"(%.3g)",get_rtime()-t_0); get_eg(&tnf1); add_eg(&eg_nf,&tnf0,&tnf1); } else @@ -1613,6 +1808,7 @@ skip_nf: get_eg(&tpz0); prim_part(nf,0,&h); get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1); + add_denomlist(BDY(h)->c); get_eg(&tnp0); if ( Demand && skip_nf_flag ) nh = newps_nosave(h,m,subst); @@ -1625,8 +1821,12 @@ skip_nf: g = updbase(g,nh); gall = append_one(gall,nh); if ( !dp_fcoeffs && ShowMag ) { - for ( mag = 0, mp = BDY(h); mp; mp = NEXT(mp) ) - mag += p_mag((P)mp->c); + for ( mag = 0, mag0 = 0, mp = BDY(h); mp; mp = NEXT(mp) ) { + magt = p_mag((P)mp->c); + mag0 = MAX(mag0,magt); + mag += magt; + } + Max_coef = MAX(Max_coef,mag0); Max_mag = MAX(Max_mag,mag); } if ( DP_Print ) { @@ -1638,7 +1838,7 @@ skip_nf: l->dp1,l->dp2,length(g),length(gall),DPPlength(d), pss[nh]); if ( ShowMag ) - fprintf(asir_out,",mag=%d",mag); + fprintf(asir_out,",mag=(%d,%d)",mag,mag0); fprintf(asir_out,"\n"); fflush(asir_out); } else if ( DP_PrintShort ) { fprintf(asir_out,"+"); fflush(asir_out); @@ -1733,7 +1933,7 @@ DP_pairs updpairs( DP_pairs d, NODE /* of index */ g, if ( !NoCriB && d ) { dl = DPPlength(d); d = criterion_B( d, t ); - dl -= DPPlength(d); NBP += dl; + dl -= DPPlength(d); N_BP += dl; } d1 = newpairs( g, t ); if ( NEXT(d1) ) { @@ -1779,7 +1979,7 @@ DP_pairs newpairs( NODE /* of index */ g, int t ) p->lcm = lcm_of_DL(CNVars, dl = psh[dp], tdl, (DL)0 ); #if 0 if ( do_weyl ) - p->sugar = dl_weight(p->lcm); + p->sugar = dl_weyl_weight(p->lcm); else #endif p->sugar = (ts > (s = pss[dp] - dl->td) ? ts : s) + p->lcm->td; @@ -1895,7 +2095,7 @@ DP_pairs criterion_F( DP_pairs d1 ) for ( head = last = 0, p = d1; NEXT(p); ) { s = (r = w = collect_pairs_of_hdlcm( p, &rest ))->sugar; while ( w = NEXT(w) ) - if ( criterion_2( w->dp1, w->dp2 ) ) { + if ( !do_weyl && criterion_2( w->dp1, w->dp2 ) ) { r = w; break; } else if ( w->sugar < s ) s = (r = w)->sugar; @@ -1946,7 +2146,8 @@ int dl_redble(DL dl1,DL dl2) return 1; } -int dl_weight(DL dl) +#if 0 +int dl_weyl_weight(DL dl) { int n,w,i; @@ -1955,6 +2156,7 @@ int dl_weight(DL dl) w += (-dl->d[i]+dl->d[n+i]); return w; } +#endif int gbcheck(NODE f) { @@ -1994,7 +2196,7 @@ int gbcheck(NODE f) return 1; } -void gbcheck_list(NODE f,VECT *gp,LIST *pp) +void gbcheck_list(NODE f,int n,VECT *gp,LIST *pp) { int i; NODE r,g,gall,u,u0,t; @@ -2003,6 +2205,8 @@ void gbcheck_list(NODE f,VECT *gp,LIST *pp) DP_pairs d,l; Q q1,q2; + /* we need the following settings */ + NVars = CNVars = n; setup_arrays(f,0,&r); for ( gall = g = 0, d = 0; r; r = NEXT(r) ) { i = (int)BDY(r); @@ -2057,6 +2261,7 @@ void dp_set_flag(Obj name,Obj value) { char *n; int v; + Q ratio; if ( OID(name) != O_STR ) return; @@ -2067,6 +2272,16 @@ void dp_set_flag(Obj name,Obj value) if ( !strcmp(n,"Dist") ) { Dist = (LIST)value; return; } + if ( !strcmp(n,"Content") ) { + ratio = (Q)value; + if ( ratio ) { + DP_Multiple = BD(NM(ratio))[0]; + Denominator = INT(ratio) ? 1 : BD(DN(ratio))[0]; + } else { + DP_Multiple = 0; + Denominator = 1; + } + } if ( value && OID(value) != O_N ) return; v = QTOS((Q)value); @@ -2106,16 +2321,26 @@ void dp_set_flag(Obj name,Obj value) GenTrace = v; else if ( !strcmp(n,"OXCheck") ) OXCheck = v; + else if ( !strcmp(n,"GenSyz") ) + GenSyz = v; } void dp_make_flaglist(LIST *list) { - Q v; + Q v,nm,dn; STRING name,path; NODE n,n1; +#if 0 STOQ(DP_Multiple,v); MKNODE(n,v,0); MKSTR(name,"DP_Multiple"); MKNODE(n1,name,n); n = n1; STOQ(Denominator,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Denominator"); MKNODE(n1,name,n); n = n1; +#else + if ( DP_Multiple ) { + STOQ(DP_Multiple,nm); STOQ(Denominator,dn); divq(nm,dn,&v); + } else + v = 0; + MKNODE(n,v,0); MKSTR(name,"Content"); MKNODE(n1,name,n); n = n1; +#endif MKNODE(n1,Dist,n); n = n1; MKSTR(name,"Dist"); MKNODE(n1,name,n); n = n1; STOQ(Reverse,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Reverse"); MKNODE(n1,name,n); n = n1; STOQ(Stat,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Stat"); MKNODE(n1,name,n); n = n1; @@ -2124,6 +2349,7 @@ void dp_make_flaglist(LIST *list) STOQ(DP_NFStat,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NFStat"); MKNODE(n1,name,n); n = n1; STOQ(OXCheck,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"OXCheck"); MKNODE(n1,name,n); n = n1; STOQ(GenTrace,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"GenTrace"); MKNODE(n1,name,n); n = n1; + STOQ(GenSyz,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"GenSyz"); MKNODE(n1,name,n); n = n1; STOQ(PtozpRA,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"PtozpRA"); MKNODE(n1,name,n); n = n1; STOQ(ShowMag,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"ShowMag"); MKNODE(n1,name,n); n = n1; STOQ(Top,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Top"); MKNODE(n1,name,n); n = n1; @@ -2168,7 +2394,11 @@ void dp_load(int index,DP *p) sprintf(path,"%s%c%d",Demand,DELIM,index); if ( !(fp = fopen(path,"rb") ) ) error("dp_load : cannot open a file"); - skipvl(fp); loadobj(fp,(Obj *)p); fclose(fp); + if ( PCoeffs ) + loadvl(fp); + else + skipvl(fp); + loadobj(fp,(Obj *)p); fclose(fp); } } @@ -2181,7 +2411,11 @@ int dp_load_t(int index,DP *p) if ( !(fp = fopen(path,"rb") ) ) return 0; else { - skipvl(fp); loadobj(fp,(Obj *)p); fclose(fp); return 1; + if ( PCoeffs ) + loadvl(fp); + else + skipvl(fp); + loadobj(fp,(Obj *)p); fclose(fp); return 1; } } @@ -2189,7 +2423,7 @@ void init_stat() { init_eg(&eg_nf); init_eg(&eg_nfm); init_eg(&eg_znfm); init_eg(&eg_pz); init_eg(&eg_np); init_eg(&eg_ra); init_eg(&eg_mc); init_eg(&eg_gc); - ZR = NZR = TP = NMP = NBP = NFP = NDP = 0; + ZR = NZR = TP = NMP = N_BP = NFP = NDP = 0; } void print_stat() { @@ -2198,7 +2432,7 @@ void print_stat() { print_eg("NF",&eg_nf); print_eg("NFM",&eg_nfm); print_eg("ZNFM",&eg_znfm); print_eg("PZ",&eg_pz); print_eg("NP",&eg_np); print_eg("RA",&eg_ra); print_eg("MC",&eg_mc); print_eg("GC",&eg_gc); - fprintf(asir_out,"T=%d,B=%d M=%d F=%d D=%d ZR=%d NZR=%d\n",TP,NBP,NMP,NFP,NDP,ZR,NZR); + fprintf(asir_out,"T=%d,B=%d M=%d F=%d D=%d ZR=%d NZR=%d\n",TP,N_BP,NMP,NFP,NDP,ZR,NZR); } /*