=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2018/engine/nd.c,v retrieving revision 1.45 retrieving revision 1.53 diff -u -p -r1.45 -r1.53 --- OpenXM_contrib2/asir2018/engine/nd.c 2021/01/11 08:37:44 1.45 +++ OpenXM_contrib2/asir2018/engine/nd.c 2021/03/12 01:18:33 1.53 @@ -1,10 +1,11 @@ -/* $OpenXM: OpenXM_contrib2/asir2018/engine/nd.c,v 1.44 2020/12/15 07:40:09 noro Exp $ */ +/* $OpenXM: OpenXM_contrib2/asir2018/engine/nd.c,v 1.52 2021/03/09 08:48:50 noro Exp $ */ #include "nd.h" void print_siglist(NODE l); -int Nnd_add,Nf4_red,NcriB,NcriMF,Ncri2,Npairs; +NODE nd_hpdata; +int Nnd_add,Nf4_red,NcriB,NcriMF,Ncri2,Npairs,Nnewpair; struct oEGT eg_search,f4_symb,f4_conv,f4_elim1,f4_elim2; int diag_period = 6; @@ -69,7 +70,7 @@ static NODE nd_tracelist; static NODE nd_alltracelist; static int nd_gentrace,nd_gensyz,nd_nora,nd_newelim,nd_intersect,nd_lf,nd_norb; static int nd_f4_td,nd_sba_f4step,nd_sba_pot,nd_sba_largelcm,nd_sba_dontsort,nd_sba_redundant_check; -static int nd_top,nd_sba_syz,nd_sba_fixord,nd_sba_grevlexgb; +static int nd_top,nd_sba_syz,nd_sba_inputisgb; static int *nd_gbblock; static NODE nd_nzlist,nd_check_splist; static int nd_splist; @@ -78,6 +79,21 @@ static int nd_f4red,nd_rank0,nd_last_nonzero; static DL *nd_sba_hm; static NODE *nd_sba_pos; +struct comp_sig_spec { + int n; + // current_i <-> oldv[i] + int *oldv; + int *weight; + struct order_pair *order_pair; + int block_length; + int **matrix; + int row; + int (*cmpdl)(int n,DL d1,DL d2); +}; + +struct comp_sig_spec *nd_sba_modord; + +DL ndltodl(int n,UINT *ndl); NumberField get_numberfield(); UINT *nd_det_compute_bound(NDV **dm,int n,int j); void nd_det_reconstruct(NDV **dm,int n,int j,NDV d); @@ -95,7 +111,7 @@ P ndc_div(int mod,union oNDC a,union oNDC b); P ndctop(int mod,union oNDC c); void finalize_tracelist(int i,P cont); void conv_ilist(int demand,int trace,NODE g,int **indp); -void parse_nd_option(NODE opt); +void parse_nd_option(VL vl,NODE opt); void dltondl(int n,DL dl,UINT *r); DP ndvtodp(int mod,NDV p); DP ndtodp(int mod,ND p); @@ -2453,6 +2469,130 @@ LIST compute_splist() return l0; } +typedef struct oHPDATA { + P hn; // HP(t)=hn(t)/(1-t)^n + int len; + P *head; // hp(i)=head[i] (i=0,...,len-1) + P hp; // dim Hm(i)=hp(i) (i >= len) + VECT x; // BDY(x)[i] = <<0,...,1,...,0>> + P *plist; // plist[i]=(1-t)^i +} *HPDATA; + +void make_reduced(VECT b,int nv); +void mhp_rec(VECT b,VECT x,P t,P *r); +P mhp_ctop(P *r,P *plist,int n); +void mhp_to_hf(VL vl,P hp,int n,P *plist,VECT *head,P *hf); +DL monomial_colon(DL a,DL b,int n); +LIST dp_monomial_hilbert_poincare(VECT b,VECT x,P *plist); + +int hpvalue(HPDATA data,int d) +{ + P *head; + int len; + P hp,val; + Z dz; + + head = data->head; + len = data->len; + hp = data->hp; + if ( d < len ) + return ZTOS((Z)head[d]); + else { + STOZ(d,dz); + substp(CO,hp,hp->v,(P)dz,&val); + return ZTOS((Z)val); + } +} + +void setup_hpdata(HPDATA final,HPDATA current) +{ + int n,i; + P *r; + DL *p; + P tv; + VECT b,x,head; + DL dl; + + n = nd_nvar; + final->hn = (P)ARG0(nd_hpdata); + head = (VECT)ARG2(nd_hpdata); + final->len = head->len; + final->head = (P *)BDY(head); + final->hp = (P)ARG3(nd_hpdata); + final->plist = (P *)BDY((VECT)ARG4(nd_hpdata)); + MKVECT(x,n); + for ( i = 0; i < n; i++ ) { + NEWDL(dl,n); dl->d[i] = 1; dl->td = 1; BDY(x)[i] = dl; + } + final->x = x; + + r = (P *)CALLOC(n+1,sizeof(P)); + MKVECT(b,nd_psn); p = (DL *)BDY(b); + for ( i = 0; i < nd_psn; i++ ) { + p[i] = ndltodl(n,nd_psh[i]->dl); + } + make_reduced(b,n); + makevar("t",&tv); + mhp_rec(b,x,tv,r); + current->hn = mhp_ctop(r,final->plist,n); + mhp_to_hf(CO,current->hn,n,final->plist,&head,¤t->hp); + current->head = (P *)BDY(head); + current->len = head->len; + current->x = x; + current->plist = final->plist; +} + +void update_hpdata(HPDATA current,int nh,int do_hf) +{ + NODE data1,nd,t; + DL new,dl; + int len,i,n; + Z dz; + DL *p; + VECT b,head; + P tv,td,s,hn,hpoly; + LIST list1; + + n = nd_nvar; + new = ndltodl(n,nd_psh[nh]->dl); + MKVECT(b,nh); p = (DL *)BDY(b); + for ( i = 0; i < nh; i++ ) { + p[i] = monomial_colon(ndltodl(n,nd_psh[i]->dl),new,n); + } + // compute HP(I:new) + list1 = dp_monomial_hilbert_poincare(b,current->x,current->plist); + data1 = BDY((LIST)list1); + // HP(I+) = H(I)-t^d*H(I:new), d=tdeg(new) + makevar("t",&tv); UTOZ(new->td,dz); + pwrp(CO,tv,dz,&td); + mulp(CO,(P)ARG0(data1),td,&s); + subp(CO,current->hn,s,&hn); + current->hn = hn; + if ( do_hf ) { + mhp_to_hf(CO,hn,n,current->plist,&head,&hpoly); + current->head = (P *)BDY(head); + current->len = head->len; + current->hp = hpoly; + } +} + +ND_pairs nd_remove_same_sugar( ND_pairs d, int sugar) +{ + struct oND_pairs root; + ND_pairs prev,cur; + + root.next = d; + prev = &root; cur = d; + while ( cur ) { + if ( SG(cur) == sugar ) + prev->next = cur->next; + else + prev = cur; + cur = cur->next; + } + return root.next; +} + /* return value = 0 => input is not a GB */ NODE nd_gb(int m,int ishomo,int checkonly,int gensyz,int **indp) @@ -2469,6 +2609,9 @@ NODE nd_gb(int m,int ishomo,int checkonly,int gensyz,i int Nnfnz = 0,Nnfz = 0; P cont; LIST list; + struct oHPDATA current_hpdata,final_hpdata; + int final_hpvalue; + struct oEGT eg1,eg2,eg_update; init_eg(&eg_update); @@ -2479,6 +2622,10 @@ init_eg(&eg_update); g = update_base(g,i); } sugar = 0; + if ( nd_hpdata ) { + if ( DP_Print ) fprintf(asir_out,"Hilbert driven algorithm.\n"); + setup_hpdata(&final_hpdata,¤t_hpdata); + } while ( d ) { again: l = nd_minp(d,&d); @@ -2495,6 +2642,18 @@ again: } sugar = SG(l); if ( DP_Print ) fprintf(asir_out,"%d",sugar); + if ( nd_hpdata ) { + if ( !compp(CO,final_hpdata.hn,current_hpdata.hn) ) + break; + else { + final_hpvalue = hpvalue(&final_hpdata,sugar); + if ( final_hpvalue == hpvalue(¤t_hpdata,sugar) ) { +// if ( DP_Print ) fprintf(asir_out,"done.\n",sugar); + d = nd_remove_same_sugar(d,sugar); + continue; + } + } + } } stat = nd_sp(m,0,l,&h); if ( !stat ) { @@ -2552,6 +2711,13 @@ get_eg(&eg1); get_eg(&eg2); add_eg(&eg_update,&eg1,&eg2); g = update_base(g,nh); FREENDP(l); + if ( nd_hpdata ) { + update_hpdata(¤t_hpdata,nh,1); + if ( final_hpvalue == hpvalue(¤t_hpdata,sugar) ) { +// if ( DP_Print ) fprintf(asir_out,"sugar=%d done.\n",sugar); + d = nd_remove_same_sugar(d,sugar); + } + } } else { Nnfz++; if ( nd_gentrace && gensyz ) { @@ -2578,7 +2744,9 @@ get_eg(&eg2); add_eg(&eg_update,&eg1,&eg2); } ND_pairs update_pairs_s(ND_pairs d,int t,NODE *syz); +int update_pairs_array_s(ND_pairs *d,int t,NODE *syz); ND_pairs nd_newpairs_s(int t ,NODE *syz); +ND_pairs *nd_newpairs_array_s(int t ,NODE *syz); int nd_nf_pbucket_s(int mod,ND g,NDV *ps,int full,ND *nf); int nd_nf_s(int mod,ND d,ND g,NDV *ps,int full,ND *nf); @@ -2586,7 +2754,38 @@ int nd_nf_s(int mod,ND d,ND g,NDV *ps,int full,ND *nf) void _copydl(int n,DL d1,DL d2); void _subfromdl(int n,DL d1,DL d2); extern int (*cmpdl)(int n,DL d1,DL d2); +int _dl_redble_ext(DL,DL,DL,int); +int primitive_irred(ND p,SIG sig) +{ + static int wpd=0,dlen=0; + static DL dquo,squo; + static UINT *quo; + int i; + + if ( dlen < nd_nvar ) { + NEWDL(dquo,nd_nvar); + NEWDL(squo,nd_nvar); + dlen = nd_nvar; + } + if ( wpd != nd_wpd ) { + wpd = nd_wpd; + quo = (UINT *)MALLOC(wpd*sizeof(UINT)); + } + for ( i = 0; i < nd_psn; i++ ) { + if ( sig->pos == nd_psh[i]->sig->pos && + _dl_redble_ext(DL(nd_psh[i]->sig),DL(sig),squo,nd_nvar) ) + if ( ndl_reducible(HDL(p),DL(nd_psh[i])) ) { + if ( DP_Print ) fprintf(asir_out,"D"); + ndl_sub(HDL(p),DL(nd_psh[i]),quo); + _ndltodl(quo,dquo); + if ( _eqdl(nd_nvar,squo,dquo) ) + return 0; + } + } + return 1; +} + NODE insert_sig(NODE l,SIG s) { int pos; @@ -2639,8 +2838,6 @@ ND_pairs remove_spair_s(ND_pairs d,SIG sig) return (ND_pairs)root.next; } -int _dl_redble_ext(DL,DL,DL,int); - int small_lcm(ND_pairs l) { SIG sig; @@ -2707,12 +2904,14 @@ ND_pairs find_smallest_lcm(ND_pairs l) _addtodl(nd_nvar,quo,mul); if ( (*cmpdl)(nd_nvar,minlm,mul) > 0 ) { minindex = i; + break; _copydl(nd_nvar,mul,minlm); } } } // l->lcm is minimal; return l itself if ( minindex < 0 ) return l; + else return 0; for ( i = 0; i < nd_psn; i++ ) { if ( i == minindex ) continue; _ndltodl(DL(nd_psh[i]),mul); @@ -2779,7 +2978,7 @@ SIG trivial_sig(int i,int j) if ( nvar != nd_nvar ) { nvar = nd_nvar; NEWDL(lcm,nvar); NEWDL(sigi.dl,nvar); NEWDL(sigj.dl,nvar); } - if ( nd_sba_grevlexgb != 0 ) { + if ( nd_sba_inputisgb != 0 ) { lcm_of_DL(nd_nvar,nd_sba_hm[i],nd_sba_hm[j],lcm); sigi.pos = i; _subdl(nd_nvar,lcm,nd_sba_hm[i],sigi.dl); sigj.pos = j; _subdl(nd_nvar,lcm,nd_sba_hm[j],sigj.dl); @@ -2794,11 +2993,35 @@ SIG trivial_sig(int i,int j) return sig; } +int nd_minsig(ND_pairs *d) +{ + int min,i,ret; + + min = -1; + for ( i = 0; i < nd_nbase; i++ ) { + if ( d[i] != 0 ) { + if ( min < 0 ) min = i; + else { + ret = comp_sig(d[i]->sig,d[min]->sig); + if ( ret < 0 ) min = i; + } + } + } + return min; +} + +int dlength(ND_pairs d) +{ + int i; + for ( i = 0; d; d = d->next, i++ ); + return i; +} + NODE nd_sba_buch(int m,int ishomo,int **indp,NODE *syzp) { int i,j,nh,sugar,stat,pos; NODE r,t,g; - ND_pairs d; + ND_pairs *d; ND_pairs l,l1; ND h,nf,s,head,nf1; NDV nfv; @@ -2808,68 +3031,68 @@ NODE nd_sba_buch(int m,int ishomo,int **indp,NODE *syz LIST list; SIG sig; NODE *syzlist; - int ngen; + int ngen,ind; int Nnominimal,Nredundant; DL lcm,quo,mul; - struct oEGT eg1,eg2,eg_update,eg_remove,eg_large,eg_nf,eg_nfzero; - int Nnfs=0,Nnfz=0,Nnfnz=0; + struct oHPDATA final_hpdata,current_hpdata; + struct oEGT eg1,eg2,eg3,eg4,eg_update,eg_remove,eg_large,eg_nf,eg_nfzero; + struct oEGT eg_minsig,eg_smallest,eg_removecont,eg_hpdata,eg_updatepairs,eg_sbabuch,eg_sp; + int Nnfs=0,Nnfz=0,Nnfnz=0,dlen,nsyz; init_eg(&eg_remove); syzlist = (NODE *)MALLOC(nd_psn*sizeof(NODE)); + d = (ND_pairs *)MALLOC(nd_psn*sizeof(ND_pairs)); + nd_nbase = nd_psn; Nsyz = 0; Nnd_add = 0; Nnominimal = 0; Nredundant = 0; - d = 0; ngen = nd_psn; - for ( i = 0; i < nd_psn; i++ ) - for ( j = i+1; j < nd_psn; j++ ) { - sig = trivial_sig(i,j); - syzlist[sig->pos] = insert_sig(syzlist[sig->pos],sig); + if ( !do_weyl ) { + for ( i = 0; i < nd_psn; i++ ) + for ( j = i+1; j < nd_psn; j++ ) { + sig = trivial_sig(i,j); + syzlist[sig->pos] = insert_sig(syzlist[sig->pos],sig); + } } + dlen = 0; for ( i = 0; i < nd_psn; i++ ) { - d = update_pairs_s(d,i,syzlist); + dlen += update_pairs_array_s(d,i,syzlist); } sugar = 0; pos = 0; + if ( nd_hpdata ) { + setup_hpdata(&final_hpdata,¤t_hpdata); + } NEWDL(lcm,nd_nvar); NEWDL(quo,nd_nvar); NEWDL(mul,nd_nvar); +init_eg(&eg_sp); init_eg(&eg_create); init_eg(&eg_merge); +init_eg(&eg_minsig); +init_eg(&eg_smallest); init_eg(&eg_large); init_eg(&eg_nf); init_eg(&eg_nfzero); - while ( d ) { -again: - if ( DP_Print ) { - int len; - ND_pairs td; - for ( td = d, len=0; td; td = td->next, len++) - ; - if ( !(len%100) ) fprintf(asir_out,"(%d)",len); - } - l = d; d = d->next; -#if 0 - if ( small_lcm(l) ) { - if ( DP_Print ) fprintf(asir_out,"M"); - Nnominimal++; - continue; - } - if ( SG(l) != sugar ) { - sugar = SG(l); - if ( DP_Print ) fprintf(asir_out,"%d",sugar); - } - sig = l->sig; - if ( DP_Print && nd_sba_pot ) { - if ( sig->pos != pos ) { - fprintf(asir_out,"[%d]",sig->pos); - pos = sig->pos; - } - } - stat = nd_sp(m,0,l,&h); -#else +init_eg(&eg_removecont); +init_eg(&eg_updatepairs); +init_eg(&eg_hpdata); +init_eg(&eg_sbabuch); +get_eg(&eg3); + while ( 1 ) { + if ( DP_Print && !nd_hpdata && dlen%1000 == 0 ) fprintf(asir_out,"(%d)",dlen); +again : +get_eg(&eg1); + ind = nd_minsig(d); +get_eg(&eg2); add_eg(&eg_minsig,&eg1,&eg2); + if ( ind < 0 ) break; + l = d[ind]; +// printf("(%d,%d)",l->i1,l->i2); print_sig(l->sig); printf("\n"); +get_eg(&eg1); l1 = find_smallest_lcm(l); +get_eg(&eg2); add_eg(&eg_smallest,&eg1,&eg2); if ( l1 == 0 ) { - if ( DP_Print ) fprintf(asir_out,"M"); + d[ind] = d[ind]->next; dlen--; +// if ( DP_Print && !nd_hpdata ) fprintf(asir_out,"M"); Nnominimal++; continue; } @@ -2884,11 +3107,11 @@ again: pos = sig->pos; } } +get_eg(&eg1); stat = nd_sp(m,0,l1,&h); -#endif +get_eg(&eg2); add_eg(&eg_sp,&eg1,&eg2); if ( !stat ) { - NEXT(l) = d; d = l; - d = nd_reconstruct(0,d); + nd_reconstruct_s(0,d); goto again; } get_eg(&eg1); @@ -2899,14 +3122,14 @@ get_eg(&eg1); #endif get_eg(&eg2); if ( !stat ) { - NEXT(l) = d; d = l; - d = nd_reconstruct(0,d); + nd_reconstruct_s(0,d); goto again; } else if ( stat == -1 ) { + d[ind] = d[ind]->next; dlen--; Nnfs++; if ( DP_Print ) { printf("S"); fflush(stdout); } - FREENDP(l); } else if ( nf ) { + d[ind] = d[ind]->next; dlen--; Nnfnz++; if ( DP_Print ) { if ( nd_sba_redundant_check ) { @@ -2921,41 +3144,70 @@ get_eg(&eg2); } add_eg(&eg_nf,&eg1,&eg2); hc = HCU(nf); + get_eg(&eg1); nd_removecont(m,nf); + get_eg(&eg2); add_eg(&eg_removecont,&eg1,&eg2); nfv = ndtondv(m,nf); nd_free(nf); nh = ndv_newps(m,nfv,0); - - d = update_pairs_s(d,nh,syzlist); + + get_eg(&eg1); + dlen += update_pairs_array_s(d,nh,syzlist); + get_eg(&eg2); add_eg(&eg_updatepairs,&eg1,&eg2); nd_sba_pos[sig->pos] = append_one(nd_sba_pos[sig->pos],nh); - FREENDP(l); + if ( nd_hpdata ) { + get_eg(&eg1); + update_hpdata(¤t_hpdata,nh,0); + get_eg(&eg2); add_eg(&eg_hpdata,&eg1,&eg2); + if ( !compp(CO,final_hpdata.hn,current_hpdata.hn) ) { + if ( DP_Print ) { printf("\nWe found a gb.\n"); } + break; + } + } } else { + d[ind] = d[ind]->next; dlen--; Nnfz++; add_eg(&eg_nfzero,&eg1,&eg2); // syzygy get_eg(&eg1); - d = remove_spair_s(d,sig); + nsyz = Nsyz; + d[sig->pos] = remove_spair_s(d[sig->pos],sig); + dlen -= Nsyz-nsyz; get_eg(&eg2); add_eg(&eg_remove,&eg1,&eg2); syzlist[sig->pos] = insert_sig(syzlist[sig->pos],sig); if ( DP_Print ) { printf("."); fflush(stdout); } - FREENDP(l); } } + get_eg(&eg4); add_eg(&eg_sbabuch,&eg3,&eg4); g = conv_ilist_s(nd_demand,0,indp); if ( DP_Print ) { - printf("\nnd_sba done. nd_add=%d,Nsyz=%d,Nsamesig=%d,Nnominimal=%d\n",Nnd_add,Nsyz,Nsamesig,Nnominimal); + printf("\ndlen=%d,nd_sba done. nd_add=%d,Nsyz=%d,Nsamesig=%d,Nnominimal=%d\n",dlen,Nnd_add,Nsyz,Nsamesig,Nnominimal); printf("Nnfnz=%d,Nnfz=%d,Nnfsingular=%d\n",Nnfnz,Nnfz,Nnfs); - fflush(stdout); + fflush(stdout); if ( nd_sba_redundant_check ) printf("Nredundant=%d\n",Nredundant); - fflush(stdout); + fflush(stdout); + print_eg("sp",&eg_sp); print_eg("create",&eg_create); print_eg("merge",&eg_merge); + print_eg("minsig",&eg_minsig); + print_eg("smallest",&eg_smallest); print_eg("remove",&eg_remove); + printf("\n"); print_eg("nf",&eg_nf); print_eg("nfzero",&eg_nfzero); + print_eg("removecont",&eg_removecont); + print_eg("updatepairs",&eg_updatepairs); + print_eg("hpdata",&eg_hpdata); + print_eg("total",&eg_sbabuch); printf("\n"); } if ( nd_sba_syz ) { + print_eg("remove",&eg_remove); + print_eg("nf",&eg_nf); + print_eg("nfzero",&eg_nfzero); + printf("\n"); + } + if ( nd_sba_syz ) { NODE hsyz,tsyz,prev; hsyz = 0; @@ -3136,6 +3388,8 @@ NODE nd_gb_trace(int m,int ishomo,int **indp) int diag_count = 0; P cont; LIST list; + struct oHPDATA current_hpdata,final_hpdata; + int final_hpvalue; init_eg(&eg_monic); init_eg(&eg_invdalg); @@ -3146,6 +3400,11 @@ NODE nd_gb_trace(int m,int ishomo,int **indp) g = update_base(g,i); } sugar = 0; + if ( nd_hpdata ) { + if ( DP_Print ) fprintf(asir_out,"Hilbert driven algorithm.\n"); + setup_hpdata(&final_hpdata,¤t_hpdata); + } + while ( d ) { again: l = nd_minp(d,&d); @@ -3166,6 +3425,18 @@ again: #endif sugar = SG(l); if ( DP_Print ) fprintf(asir_out,"%d",sugar); + if ( nd_hpdata ) { + if ( !compp(CO,final_hpdata.hn,current_hpdata.hn) ) + break; + else { + final_hpvalue = hpvalue(&final_hpdata,sugar); + if ( final_hpvalue == hpvalue(¤t_hpdata,sugar) ) { +// if ( DP_Print ) fprintf(asir_out,"sugar=%d done.\n",sugar); + d = nd_remove_same_sugar(d,sugar); + continue; + } + } + } } stat = nd_sp(m,0,l,&h); if ( !stat ) { @@ -3239,6 +3510,13 @@ again: } d = update_pairs(d,g,nh,0); g = update_base(g,nh); + if ( nd_hpdata ) { + update_hpdata(¤t_hpdata,nh,1); + if ( final_hpvalue == hpvalue(¤t_hpdata,sugar) ) { +// if ( DP_Print ) fprintf(asir_out,"sugar=%d done.\n",sugar); + d = nd_remove_same_sugar(d,sugar); + } + } } else { if ( DP_Print ) { printf("*"); fflush(stdout); } } @@ -3415,6 +3693,23 @@ get_eg(&eg3); add_eg(&eg_merge,&eg2,&eg3); return d; } +int update_pairs_array_s( ND_pairs *d, int t,NODE *syz) +{ + ND_pairs *d1; + struct oEGT eg1,eg2,eg3; + int i; + + if ( !t ) return 0; +get_eg(&eg1); + Nnewpair = 0; + d1 = nd_newpairs_array_s(t,syz); +get_eg(&eg2); add_eg(&eg_create,&eg1,&eg2); + for ( i = 0; i < nd_nbase; i++ ) + d[i] = merge_pairs_s(d[i],d1[i]); +get_eg(&eg3); add_eg(&eg_merge,&eg2,&eg3); + return Nnewpair; +} + ND_pairs nd_newpairs( NODE g, int t ) { NODE h; @@ -3452,12 +3747,205 @@ ND_pairs nd_newpairs( NODE g, int t ) return r0; } +int sig_cmpdl_op(int n,DL d1,DL d2) +{ + int e1,e2,i,j,l; + int *t1,*t2; + int len,head; + struct order_pair *pair; + + len = nd_sba_modord->block_length; + pair = nd_sba_modord->order_pair; + + head = 0; + for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) { + l = pair[i].length; + switch ( pair[i].order ) { + case 0: + for ( j = 0, e1 = e2 = 0; j < l; j++ ) { + e1 += t1[j]; + e2 += t2[j]; + } + if ( e1 > e2 ) + return 1; + else if ( e1 < e2 ) + return -1; + else { + for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- ); + if ( j >= 0 ) + return t1[j] < t2[j] ? 1 : -1; + } + break; + case 1: + for ( j = 0, e1 = e2 = 0; j < l; j++ ) { + e1 += t1[j]; + e2 += t2[j]; + } + if ( e1 > e2 ) + return 1; + else if ( e1 < e2 ) + return -1; + else { + for ( j = 0; j < l && t1[j] == t2[j]; j++ ); + if ( j < l ) + return t1[j] > t2[j] ? 1 : -1; + } + break; + case 2: + for ( j = 0; j < l && t1[j] == t2[j]; j++ ); + if ( j < l ) + return t1[j] > t2[j] ? 1 : -1; + break; + default: + error("sig_cmpdl_op : invalid order"); break; + } + t1 += l; t2 += l; head += l; + } + return 0; +} + +int sig_cmpdl_mat(int n,DL d1,DL d2) +{ + int *v,*t1,*t2; + int s,i,j,len; + int **matrix; + static int *w; + static int nvar = 0; + + if ( nvar != n ) { + nvar = n; w = (int *)MALLOC(n*sizeof(int)); + } + for ( i = 0, t1 = d1->d, t2 = d2->d; i < n; i++ ) + w[i] = t1[i]-t2[i]; + len = nd_sba_modord->row; + matrix = nd_sba_modord->matrix; + for ( j = 0; j < len; j++ ) { + v = matrix[j]; + for ( i = 0, s = 0; i < n; i++ ) + s += v[i]*w[i]; + if ( s > 0 ) + return 1; + else if ( s < 0 ) + return -1; + } + return 0; +} + +struct comp_sig_spec *create_comp_sig_spec(VL current_vl,VL old_vl,Obj ord,Obj weight) +{ + struct comp_sig_spec *spec; + VL ovl,vl; + V ov; + int i,j,n,nvar,s; + NODE node,t,tn; + struct order_pair *l; + MAT m; + Obj **b; + int **w; + int *a; + + spec = (struct comp_sig_spec *)MALLOC(sizeof(struct comp_sig_spec)); + for ( i = 0, vl = current_vl; vl; vl = NEXT(vl), i++ ); + spec->n = nvar = i; + if ( old_vl != 0 ) { + spec->oldv = (int *)MALLOC(nvar*sizeof(int)); + for ( i = 0, ovl = old_vl; i < nvar; ovl = NEXT(ovl), i++ ) { + ov = ovl->v; + for ( j = 0, vl = current_vl; vl; vl = NEXT(vl), j++ ) + if ( ov == vl->v ) break; + spec->oldv[i] = j; + } + } else + spec->oldv = 0; + if ( !ord || NUM(ord) ) { + switch ( ZTOS((Z)ord) ) { + case 0: + spec->cmpdl = cmpdl_revgradlex; break; + case 1: + spec->cmpdl = cmpdl_gradlex; break; + case 2: + spec->cmpdl = cmpdl_lex; break; + default: + error("create_comp_sig_spec : invalid spec"); break; + } + } else if ( OID(ord) == O_LIST ) { + node = BDY((LIST)ord); + for ( n = 0, t = node; t; t = NEXT(t), n++ ); + l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair)); + for ( i = 0, t = node, s = 0; i < n; t = NEXT(t), i++ ) { + tn = BDY((LIST)BDY(t)); l[i].order = ZTOS((Q)BDY(tn)); + tn = NEXT(tn); l[i].length = ZTOS((Q)BDY(tn)); + s += l[i].length; + } + if ( s != nvar ) + error("create_comp_sig_spec : invalid spec"); + spec->order_pair = l; + spec->block_length = n; + spec->cmpdl = sig_cmpdl_op; + } else if ( OID(ord) == O_MAT ) { + m = (MAT)ord; b = (Obj **)BDY(m); + if ( m->col != nvar ) + error("create_comp_sig_spec : invalid spec"); + w = almat(m->row,m->col); + for ( i = 0; i < m->row; i++ ) + for ( j = 0; j < m->col; j++ ) + w[i][j] = ZTOS((Q)b[i][j]); + spec->row = m->row; + spec->matrix = w; + spec->cmpdl = sig_cmpdl_mat; + } else + error("create_comp_sig_spec : invalid spec"); + if ( weight != 0 ) { + node = BDY((LIST)weight); + a = (int *)MALLOC(nvar*sizeof(int)); + for ( i = 0; i < nvar; i++, node = NEXT(node) ) + a[i] = ZTOS((Z)BDY(node)); + spec->weight = a; + } + return spec; +} + +#define SIG_MUL_WEIGHT(a,i) (weight?(a)*weight[i]:(a)) + +int comp_sig_monomial(int n,DL d1,DL d2) +{ + static DL m1,m2; + static int nvar = 0; + int *oldv,*weight; + int i,w1,w2; + + if ( nvar != n ) { + nvar = n; NEWDL(m1,nvar); NEWDL(m2,nvar); + } + if ( !nd_sba_modord ) + return (*cmpdl)(n,d1,d2); + else { + weight = nd_sba_modord->weight; + oldv = nd_sba_modord->oldv; + if ( oldv ) { + for ( i = 0; i < n; i++ ) { + m1->d[i] = d1->d[oldv[i]]; m2->d[i] = d2->d[oldv[i]]; + } + } else { + for ( i = 0; i < n; i++ ) { + m1->d[i] = d1->d[i]; m2->d[i] = d2->d[i]; + } + } + for ( i = 0, w1 = w2 = 0; i < n; i++ ) { + w1 += SIG_MUL_WEIGHT(m1->d[i],i); + w2 += SIG_MUL_WEIGHT(m2->d[i],i); + } + m1->td = w1; m2->td = w2; + return (*nd_sba_modord->cmpdl)(n,m1,m2); + } +} + int comp_sig(SIG s1,SIG s2) { if ( nd_sba_pot ) { if ( s1->pos > s2->pos ) return 1; else if ( s1->pos < s2->pos ) return -1; - else return (*cmpdl)(nd_nvar,s1->dl,s2->dl); + else return comp_sig_monomial(nd_nvar,s1->dl,s2->dl); } else { static DL m1,m2; static int nvar = 0; @@ -3468,10 +3956,10 @@ int comp_sig(SIG s1,SIG s2) } _adddl(nd_nvar,s1->dl,nd_sba_hm[s1->pos],m1); _adddl(nd_nvar,s2->dl,nd_sba_hm[s2->pos],m2); - if ( nd_sba_fixord ) - ret = cmpdl_revgradlex(nd_nvar,m1,m2); - else + if ( !nd_sba_modord ) ret = (*cmpdl)(nd_nvar,m1,m2); + else + ret = comp_sig_monomial(nd_nvar,m1,m2); if ( ret != 0 ) return ret; else if ( s1->pos > s2->pos ) return 1; else if ( s1->pos < s2->pos ) return -1; @@ -3491,9 +3979,11 @@ int _create_spair_s(int i1,int i2,ND_pairs sp,SIG sig1 p1 = nd_psh[i1]; p2 = nd_psh[i2]; ndl_lcm(DL(p1),DL(p2),sp->lcm); +#if 0 s1 = SG(p1)-TD(DL(p1)); s2 = SG(p2)-TD(DL(p2)); SG(sp) = MAX(s1,s2) + TD(sp->lcm); +#endif if ( wpd != nd_wpd ) { wpd = nd_wpd; @@ -3519,6 +4009,11 @@ int _create_spair_s(int i1,int i2,ND_pairs sp,SIG sig1 if ( ret == 0 ) return 0; else if ( ret > 0 ) sp->sig = sig1; else sp->sig = sig2; + + s1 = DL(sig1)->td+nd_sba_hm[p1->sig->pos]->td; + s2 = DL(sig2)->td+nd_sba_hm[p2->sig->pos]->td; + SG(sp) = MAX(s1,s2); + return 1; } @@ -3559,6 +4054,7 @@ ND_pairs merge_pairs_s(ND_pairs p1,ND_pairs p2) } else if ( ret > 0 ) { r->next = q2; r = q2; q2 = q2->next; } else { + Nnewpair--; ret = DL_COMPARE(q1->lcm,q2->lcm); Nsamesig++; if ( ret < 0 ) { @@ -3581,7 +4077,7 @@ ND_pairs merge_pairs_s(ND_pairs p1,ND_pairs p2) ND_pairs insert_pair_s(ND_pairs l,ND_pairs s) { ND_pairs p,prev; - int ret; + int ret=1; for ( p = l, prev = 0; p != 0; prev = p, p = p->next ) { if ( (ret = comp_sig(s->sig,p->sig)) <= 0 ) @@ -3602,6 +4098,7 @@ ND_pairs insert_pair_s(ND_pairs l,ND_pairs s) return l; } else { // insert s between prev and p + Nnewpair++; s->next = p; if ( prev == 0 ) { return s; @@ -3663,6 +4160,44 @@ ND_pairs nd_newpairs_s(int t, NODE *syz) return r0; } +ND_pairs *nd_newpairs_array_s(int t, NODE *syz) +{ + NODE h,s; + UINT *dl; + int ts,ret,i; + ND_pairs r,r0,_sp,sp; + ND_pairs *d; + SIG spsig,tsig; + static int nvar = 0; + static SIG _sig1,_sig2; + struct oEGT eg1,eg2,eg3,eg4; + + NEWND_pairs(_sp); + if ( !_sig1 || nvar != nd_nvar ) { + nvar = nd_nvar; NEWSIG(_sig1); NEWSIG(_sig2); + } + d = (ND_pairs *)MALLOC(nd_nbase*sizeof(ND_pairs)); + Nnewpair = 0; + for ( i = 0; i < t; i++ ) { + ret = _create_spair_s(i,t,_sp,_sig1,_sig2); + if ( ret ) { + spsig = _sp->sig; + for ( s = syz[spsig->pos]; s; s = s->next ) { + tsig = (SIG)s->body; + if ( _dl_redble(DL(tsig),DL(spsig),nd_nvar) ) + break; + } + if ( s == 0 ) { + NEWND_pairs(sp); + dup_ND_pairs(sp,_sp); + d[spsig->pos] = insert_pair_s(d[spsig->pos],sp); + } else + Nsyz++; + } + } + return d; +} + /* ipair = [i1,i2],[i1,i2],... */ ND_pairs nd_ipairtospair(NODE ipair) { @@ -4033,8 +4568,8 @@ int ndv_newps(int m,NDV a,NDV aq) return nd_psn++; } -// find LM wrt grevlex -void ndv_lm_fixord(NDV p,DL d) +// find LM wrt the specified modord +void ndv_lm_modord(NDV p,DL d) { NMV m; DL tmp; @@ -4042,13 +4577,13 @@ void ndv_lm_fixord(NDV p,DL d) NEWDL(tmp,nd_nvar); m = BDY(p); len = LEN(p); - _ndltodl(DL(m),d); printdl(d); printf("->"); + _ndltodl(DL(m),d); // printdl(d); printf("->"); for ( i = 1, NMV_ADV(m); i < len; i++, NMV_ADV(m) ) { _ndltodl(DL(m),tmp); - ret = cmpdl_revgradlex(nd_nvar,tmp,d); + ret = comp_sig_monomial(nd_nvar,tmp,d); if ( ret > 0 ) _copydl(nd_nvar,tmp,d); } - printdl(d); printf("\n"); +// printdl(d); printf("\n"); } /* nd_tracelist = [[0,index,div],...,[nd_psn-1,index,div]] */ @@ -4169,8 +4704,8 @@ int ndv_setup(int mod,int trace,NODE f,int dont_sort,i if ( nd_demand ) nd_ps_trace_sym[i]->sig = sig; } NEWDL(nd_sba_hm[i],nd_nvar); - if ( nd_sba_fixord ) - ndv_lm_fixord(nd_ps[i],nd_sba_hm[i]); + if ( nd_sba_modord ) + ndv_lm_modord(nd_ps[i],nd_sba_hm[i]); else _ndltodl(DL(nd_psh[i]),nd_sba_hm[i]); } @@ -4318,7 +4853,6 @@ void nd_gr(LIST f,LIST v,int m,int homo,int retdp,int nd_module = 0; if ( !m && Demand ) nd_demand = 1; else nd_demand = 0; - parse_nd_option(current_option); if ( DP_Multiple ) nd_scale = ((double)DP_Multiple)/(double)(Denominator?Denominator:1); @@ -4326,6 +4860,7 @@ void nd_gr(LIST f,LIST v,int m,int homo,int retdp,int ndv_alloc = 0; #endif get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc); + parse_nd_option(vv,current_option); if ( m && nd_vc ) error("nd_{gr,f4} : computation over Fp(X) is unsupported. Use dp_gr_mod_main()."); for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); @@ -4550,15 +5085,15 @@ void nd_sba(LIST f,LIST v,int m,int homo,int retdp,int int *perm; EPOS oepos; int obpe,oadv,ompos,cbpe; - struct oEGT eg0,eg1,egconv; + struct oEGT eg0,eg1,egconv,egintred; nd_module = 0; nd_demand = 0; - parse_nd_option(current_option); Nsamesig = 0; if ( DP_Multiple ) nd_scale = ((double)DP_Multiple)/(double)(Denominator?Denominator:1); get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc); + parse_nd_option(vv,current_option); if ( m && nd_vc ) error("nd_sba : computation over Fp(X) is unsupported. Use dp_gr_mod_main()."); for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); @@ -4636,8 +5171,10 @@ void nd_sba(LIST f,LIST v,int m,int homo,int retdp,int nd_setup_parameters(nvar,0); } nd_demand = 0; + get_eg(&eg0); x = ndv_reducebase(x,perm); x = ndv_reduceall(m,x); + get_eg(&eg1); init_eg(&egintred); add_eg(&egintred,&eg0,&eg1); nd_setup_parameters(nd_nvar,0); get_eg(&eg0); for ( r0 = 0, t = x; t; t = NEXT(t) ) { @@ -4657,6 +5194,7 @@ void nd_sba(LIST f,LIST v,int m,int homo,int retdp,int } else MKLIST(*rp,r0); get_eg(&eg1); init_eg(&egconv); add_eg(&egconv,&eg0,&eg1); + print_eg("intred",&egintred); fprintf(asir_out,"\n"); print_eg("conv",&egconv); fprintf(asir_out,"\n"); } @@ -4676,8 +5214,8 @@ void nd_gr_postproc(LIST f,LIST v,int m,struct order_s struct order_spec *ord1; int *perm; - parse_nd_option(current_option); get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc); + parse_nd_option(vv,current_option); for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); switch ( ord->id ) { case 1: @@ -4824,8 +5362,8 @@ void nd_gr_recompute_trace(LIST f,LIST v,int m,struct int len,n,j; NDV *db,*pb; - parse_nd_option(current_option); get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc); + parse_nd_option(vv,current_option); for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); switch ( ord->id ) { case 1: @@ -4924,7 +5462,8 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,int NcriB = NcriMF = Ncri2 = 0; nd_module = 0; nd_lf = 0; - parse_nd_option(current_option); + get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc); + parse_nd_option(vv,current_option); if ( nd_lf ) { if ( f4 ) nd_f4_lf_trace(f,v,trace,homo,ord,rp); @@ -4935,7 +5474,6 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,int if ( DP_Multiple ) nd_scale = ((double)DP_Multiple)/(double)(Denominator?Denominator:1); - get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc); for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); switch ( ord->id ) { case 1: @@ -5966,6 +6504,7 @@ void nd_reconstruct_s(int trace,ND_pairs *d) SG(s) = SG(t); ndl_reconstruct(LCM(t),LCM(s),obpe,oepos); } + if ( s0 ) NEXT(s) = 0; d[i] = s0; } @@ -6681,7 +7220,7 @@ ND ptond(VL vl,VL dvl,P p) else if ( NUM(p) ) { NEWNM(m); ndl_zero(DL(m)); - if ( !INT((Q)p) ) + if ( RATN(p) && !INT((Q)p) ) error("ptond : input must be integer-coefficient"); CZ(m) = (Z)p; NEXT(m) = 0; @@ -10084,12 +10623,13 @@ NODE conv_ilist_s(int demand,int trace,int **indp) return g0; } -void parse_nd_option(NODE opt) +void parse_nd_option(VL vl,NODE opt) { NODE t,p,u; int i,s,n; char *key; Obj value; + VL oldvl; nd_gentrace = 0; nd_gensyz = 0; nd_nora = 0; nd_norb = 0; nd_gbblock = 0; nd_newelim = 0; nd_intersect = 0; nd_nzlist = 0; @@ -10097,7 +10637,8 @@ void parse_nd_option(NODE opt) nd_sugarweight = 0; nd_f4red =0; nd_rank0 = 0; nd_f4_td = 0; nd_sba_f4step = 2; nd_sba_pot = 0; nd_sba_largelcm = 0; nd_sba_dontsort = 0; nd_top = 0; nd_sba_redundant_check = 0; - nd_sba_syz = 0; nd_sba_fixord = 0; nd_sba_grevlexgb = 0; + nd_sba_syz = 0; nd_sba_modord = 0; nd_sba_inputisgb = 0; + nd_hpdata = 0; for ( t = opt; t; t = NEXT(t) ) { p = BDY((LIST)BDY(t)); @@ -10145,6 +10686,9 @@ void parse_nd_option(NODE opt) nd_splist = value?1:0; } else if ( !strcmp(key,"check_splist") ) { nd_check_splist = BDY((LIST)value); + } else if ( !strcmp(key,"hpdata") ) { + if ( value ) + nd_hpdata = BDY((LIST)value); } else if ( !strcmp(key,"sugarweight") ) { u = BDY((LIST)value); n = length(u); @@ -10163,10 +10707,19 @@ void parse_nd_option(NODE opt) nd_sba_dontsort = value?1:0; } else if ( !strcmp(key,"sba_syz") ) { nd_sba_syz = value?1:0; - } else if ( !strcmp(key,"sba_fixord") ) { - nd_sba_fixord = value?1:0; - } else if ( !strcmp(key,"sba_grevlexgb") ) { - nd_sba_grevlexgb = value?1:0; + } else if ( !strcmp(key,"sba_modord") ) { + // value=[vlist,ordspec,weight] + u = BDY((LIST)value); + pltovl((LIST)ARG0(u),&oldvl); + nd_sba_modord = create_comp_sig_spec(vl,oldvl,(Obj)ARG1(u),argc(u)==3?ARG2(u):0); + } else if ( !strcmp(key,"sba_gbinput") ) { + nd_sba_inputisgb = value?1:0; + if ( nd_sba_inputisgb != 0 ) { + // value=[vlist,ordspec,weight] + u = BDY((LIST)value); + pltovl((LIST)ARG0(u),&oldvl); + nd_sba_modord = create_comp_sig_spec(vl,oldvl,(Obj)ARG1(u),argc(u)==3?ARG2(u):0); + } } else if ( !strcmp(key,"sba_redundant_check") ) { nd_sba_redundant_check = value?1:0; } else if ( !strcmp(key,"top") ) { @@ -10378,8 +10931,8 @@ MAT nd_btog(LIST f,LIST v,int mod,struct order_spec *o if ( mod == -2 ) return nd_btog_lf(f,v,ord,tlist,rp); - parse_nd_option(current_option); get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc); + parse_nd_option(vv,current_option); for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); switch ( ord->id ) { case 1: @@ -10447,8 +11000,8 @@ MAT nd_btog_lf(LIST f,LIST v,struct order_spec *ord,LI LM lm; Z lf,inv; - parse_nd_option(current_option); get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc); + parse_nd_option(vv,current_option); for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); switch ( ord->id ) { case 1: @@ -10518,8 +11071,8 @@ VECT nd_btog_one(LIST f,LIST v,int mod,struct order_sp if ( mod == -2 ) error("nd_btog_one : not implemented yet for a large finite field"); - parse_nd_option(current_option); get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc); + parse_nd_option(vv,current_option); for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); switch ( ord->id ) { case 1: @@ -10632,8 +11185,8 @@ void nd_f4_lf_trace(LIST f,LIST v,int trace,int homo,s Q jq,bpe; nd_module = 0; - parse_nd_option(current_option); get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc); + parse_nd_option(vv,current_option); if ( nd_vc ) error("nd_f4_lf_trace : computation over a rational function field is not implemented"); for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ );