=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2018/builtin/dp.c,v retrieving revision 1.6 retrieving revision 1.32 diff -u -p -r1.6 -r1.32 --- OpenXM_contrib2/asir2018/builtin/dp.c 2019/03/18 07:00:33 1.6 +++ OpenXM_contrib2/asir2018/builtin/dp.c 2022/09/10 04:04:50 1.32 @@ -45,7 +45,7 @@ * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. * - * $OpenXM: OpenXM_contrib2/asir2018/builtin/dp.c,v 1.5 2019/03/13 08:01:05 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2018/builtin/dp.c,v 1.31 2021/12/05 22:41:03 noro Exp $ */ #include "ca.h" #include "base.h" @@ -59,9 +59,9 @@ extern struct order_spec *dp_current_spec; extern struct modorder_spec *dp_current_modspec; extern int nd_rref2; -int do_weyl; +extern int do_weyl; -void Pdp_monomial_hilbert_poincare(); +void Pdp_monomial_hilbert_poincare(),Pdp_monomial_hilbert_poincare_incremental(); void Pdp_sort(); void Pdp_mul_trunc(),Pdp_quo(); void Pdp_ord(), Pdp_ptod(), Pdp_dtop(), Phomogenize(); @@ -79,6 +79,7 @@ void Pdp_nf_mod(),Pdp_true_nf_mod(); void Pdp_criB(),Pdp_nelim(); void Pdp_minp(),Pdp_sp_mod(); void Pdp_homo(),Pdp_dehomo(); +void Pdpm_homo(),Pdpm_dehomo(),Pdpm_mod(); void Pdp_gr_mod_main(),Pdp_gr_f_main(); void Pdp_gr_main(),Pdp_gr_hm_main(),Pdp_gr_d_main(),Pdp_gr_flags(); void Pdp_interreduce(); @@ -89,8 +90,12 @@ void Pdp_vtoe(), Pdp_etov(), Pdp_dtov(), Pdp_idiv(), P void Pdp_cont(); void Pdp_gr_checklist(); void Pdp_ltod(),Pdpv_ord(),Pdpv_ht(),Pdpv_hm(),Pdpv_hc(); -void Pdpm_ltod(),Pdpm_dtol(),Pdpm_ord(),Pdpm_nf(),Pdpm_weyl_nf(),Pdpm_sp(),Pdpm_weyl_sp(); -void Pdpm_hm(),Pdpm_ht(),Pdpm_hc(); +void Pdpm_ltod(),Pdpm_dtol(),Pdpm_set_schreyer(),Pdpm_nf(),Pdpm_weyl_nf(),Pdpm_sp(),Pdpm_weyl_sp(),Pdpm_nf_and_quotient(),Pdpm_nf_and_quotient2(); +void Pdpm_schreyer_frame(),Pdpm_set_schreyer_level(); +void Pdpm_list_to_array(),Pdpm_sp_nf(),Pdpm_insert_to_zlist(); +void Pdpm_hm(),Pdpm_ht(),Pdpm_hc(),Pdpm_hp(),Pdpm_rest(),Pdpm_shift(),Pdpm_split(),Pdpm_extract(),Pdpm_sort(),Pdpm_dptodpm(),Pdpm_redble(); +void Pdpm_schreyer_base(),Pdpm_simplify_syz(),Pdpm_td(); +void Pdpm_remove_cont(); void Pdp_weyl_red(); void Pdp_weyl_sp(); @@ -108,6 +113,8 @@ void Pdp_nf_f(),Pdp_weyl_nf_f(); void Pdpm_nf_f(),Pdpm_weyl_nf_f(); void Pdp_lnf_f(); void Pnd_gr(),Pnd_gr_trace(),Pnd_f4(),Pnd_f4_trace(); +void Pnd_sba(),Pnd_sba_f4(); +void Pnd_weyl_sba(); void Pnd_gr_postproc(), Pnd_weyl_gr_postproc(); void Pnd_gr_recompute_trace(), Pnd_btog(); void Pnd_weyl_gr(),Pnd_weyl_gr_trace(); @@ -126,6 +133,7 @@ void Pdp_rref2(),Psumi_updatepairs(),Psumi_symbolic(); LIST dp_initial_term(); LIST dp_order(); +int peek_option(NODE opt,char *find,Obj *ret); void parse_gr_option(LIST f,NODE opt,LIST *v,Num *homo, int *modular,struct order_spec **ord); NODE dp_inv_or_split(NODE gb,DP f,struct order_spec *spec, DP *inv); @@ -139,6 +147,7 @@ struct ftab dp_tab[] = { {"dp_prim",Pdp_prim,1}, {"dp_red_coef",Pdp_red_coef,2}, {"dp_cont",Pdp_cont,1}, + {"dpm_remove_cont",Pdpm_remove_cont,1}, /* polynomial ring */ /* special operations */ @@ -157,9 +166,11 @@ struct ftab dp_tab[] = { {"dp_nf",Pdp_nf,4}, {"dp_nf_mod",Pdp_nf_mod,5}, {"dp_nf_f",Pdp_nf_f,4}, - {"dpm_nf_f",Pdpm_nf_f,4}, - {"dpm_weyl_nf_f",Pdpm_weyl_nf_f,4}, - {"dpm_nf",Pdpm_nf,4}, + {"dpm_nf_and_quotient",Pdpm_nf_and_quotient,-3}, + {"dpm_nf_and_quotient2",Pdpm_nf_and_quotient2,-3}, + {"dpm_nf_f",Pdpm_nf_f,-4}, + {"dpm_weyl_nf_f",Pdpm_weyl_nf_f,-4}, + {"dpm_nf",Pdpm_nf,-4}, {"dpm_sp",Pdpm_sp,2}, {"dpm_weyl_sp",Pdpm_weyl_sp,2}, @@ -186,6 +197,9 @@ struct ftab dp_tab[] = { {"dp_gr_checklist",Pdp_gr_checklist,2}, {"nd_f4",Pnd_f4,-4}, {"nd_gr",Pnd_gr,-4}, + {"nd_sba",Pnd_sba,-4}, + {"nd_weyl_sba",Pnd_weyl_sba,-4}, + {"nd_sba_f4",Pnd_sba_f4,-4}, {"nd_gr_trace",Pnd_gr_trace,-5}, {"nd_f4_trace",Pnd_f4_trace,-5}, {"nd_gr_postproc",Pnd_gr_postproc,5}, @@ -215,7 +229,7 @@ struct ftab dp_tab[] = { /* normal form */ {"dp_weyl_nf",Pdp_weyl_nf,4}, - {"dpm_weyl_nf",Pdpm_weyl_nf,4}, + {"dpm_weyl_nf",Pdpm_weyl_nf,-4}, {"dp_weyl_nf_mod",Pdp_weyl_nf_mod,5}, {"dp_weyl_nf_f",Pdp_weyl_nf_f,4}, @@ -236,6 +250,7 @@ struct ftab dp_tab[] = { /* Hilbert function */ {"dp_monomial_hilbert_poincare",Pdp_monomial_hilbert_poincare,2}, + {"dp_monomial_hilbert_poincare_incremental",Pdp_monomial_hilbert_poincare_incremental,3}, /* misc */ {"dp_inv_or_split",Pdp_inv_or_split,3}, @@ -252,7 +267,9 @@ struct ftab dp_supp_tab[] = { /* setting flags */ {"dp_sort",Pdp_sort,1}, {"dp_ord",Pdp_ord,-1}, - {"dpm_ord",Pdpm_ord,-1}, + {"dpm_set_schreyer",Pdpm_set_schreyer,-1}, + {"dpm_set_schreyer_level",Pdpm_set_schreyer_level,1}, + {"dpm_schreyer_frame",Pdpm_schreyer_frame,1}, {"dpv_ord",Pdpv_ord,-2}, {"dp_set_kara",Pdp_set_kara,-1}, {"dp_nelim",Pdp_nelim,-1}, @@ -274,7 +291,11 @@ struct ftab dp_supp_tab[] = { {"dp_ltod",Pdp_ltod,-2}, {"dpm_ltod",Pdpm_ltod,2}, - {"dpm_dtol",Pdpm_dtol,3}, + {"dpm_dptodpm",Pdpm_dptodpm,2}, + {"dpm_dtol",Pdpm_dtol,2}, + {"dpm_homo",Pdpm_homo,1}, + {"dpm_dehomo",Pdpm_dehomo,1}, + {"dpm_mod",Pdpm_mod,2}, /* criteria */ {"dp_cri1",Pdp_cri1,2}, @@ -293,6 +314,12 @@ struct ftab dp_supp_tab[] = { {"dpm_hm",Pdpm_hm,1}, {"dpm_ht",Pdpm_ht,1}, {"dpm_hc",Pdpm_hc,1}, + {"dpm_hp",Pdpm_hp,1}, + {"dpm_rest",Pdpm_rest,1}, + {"dpm_shift",Pdpm_shift,2}, + {"dpm_split",Pdpm_split,2}, + {"dpm_extract",Pdpm_extract,2}, + {"dpm_sort",Pdpm_sort,1}, {"dp_rest",Pdp_rest,1}, {"dp_initial_term",Pdp_initial_term,1}, {"dp_order",Pdp_order,1}, @@ -303,10 +330,12 @@ struct ftab dp_supp_tab[] = { {"dp_mag",Pdp_mag,1}, {"dp_sugar",Pdp_sugar,1}, {"dp_set_sugar",Pdp_set_sugar,2}, + {"dpm_td",Pdpm_td,1}, /* misc */ {"dp_mbase",Pdp_mbase,1}, {"dp_redble",Pdp_redble,2}, + {"dpm_redble",Pdpm_redble,2}, {"dp_sep",Pdp_sep,2}, {"dp_idiv",Pdp_idiv,2}, {"dp_tdiv",Pdp_tdiv,2}, @@ -316,6 +345,11 @@ struct ftab dp_supp_tab[] = { {"dp_compute_essential_df",Pdp_compute_essential_df,2}, {"dp_mono_raddec",Pdp_mono_raddec,2}, {"dp_mono_reduce",Pdp_mono_reduce,2}, + {"dpm_schreyer_base",Pdpm_schreyer_base,1}, + {"dpm_list_to_array",Pdpm_list_to_array,1}, + {"dpm_sp_nf",Pdpm_sp_nf,4}, + {"dpm_insert_to_zlist",Pdpm_insert_to_zlist,3}, + {"dpm_simplify_syz",Pdpm_simplify_syz,2}, {"dp_rref2",Pdp_rref2,2}, {"sumi_updatepairs",Psumi_updatepairs,3}, @@ -338,171 +372,6 @@ int comp_by_tdeg(DP *a,DP *b) else return 0; } -#if 0 -void make_reduced(VECT b) -{ - int n,i,j; - DP *p; - DP pi; - - n = b->len; - p = (DP *)BDY(b); - if ( BDY(p[0])->dl->td == 0 ) { - b->len = 1; - return; - } - for ( i = 0; i < n; i++ ) { - pi = p[i]; - if ( !pi ) continue; - for ( j = 0; j < n; j++ ) - if ( i != j && p[j] && dp_redble(p[j],pi) ) p[j] = 0; - } - for ( i = j = 0; i < n; i++ ) - if ( p[i] ) p[j++] = p[i]; - b->len = j; -} - -void make_reduced2(VECT b,int k) -{ - int n,i,j,l; - DP *p; - DP pi; - - n = b->len; - p = (DP *)BDY(b); - for ( i = l = k; i < n; i++ ) { - pi = p[i]; - for ( j = 0; j < k; j++ ) - if ( dp_redble(pi,p[j]) ) break; - if ( j == k ) - p[l++] = pi; - } - b->len = l; -} - -struct oEGT eg_comp; - -void mhp_rec(VECT b,VECT x,P t,P *r) -{ - int n,i,j,k,l,i2,y,len; - int *d; - Z mone,z; - DCP dc,dc1; - P s; - P *r2; - DP *p,*q; - DP pi,xj; - VECT c; - struct oEGT eg0,eg1; - - n = b->len; - y = x->len; - p = (DP *)BDY(b); - if ( !n ) { - r[0] = (P)ONE; - return; - } - if ( n == 1 && BDY(p[0])->dl->td == 0 ) { - return; - } - for ( i = 0; i < n; i++ ) - if ( BDY(p[i])->dl->td > 1 ) break; - if ( i == n ) { - r[n] = (P)ONE; - return; - } - get_eg(&eg0); - pi = p[i]; - d = BDY(pi)->dl->d; - for ( j = 0; j < y; j++ ) - if ( d[j] ) break; - xj = BDY(x)[j]; - - MKVECT(c,n); q = (DP *)BDY(c); - for ( i = k = l = 0; i < n; i++ ) - if ( BDY(p[i])->dl->d[j] ) - dp_subd(p[i],xj,&p[k++]); - else - q[l++] = p[i]; - for ( i = k, i2 = 0; i2 < l; i++, i2++ ) - p[i] = q[i2]; - /* b=(b[0]/xj,...,b[k-1]/xj,b[k],...b[n-1]) where - b[0],...,b[k-1] are divisible by k */ - make_reduced2(b,k); - get_eg(&eg1); add_eg(&eg_comp,&eg0,&eg1); - mhp_rec(b,x,t,r); - /* c = (b[0],...,b[l-1],xj) */ - q[l] = xj; c->len = l+1; - r2 = (P *)CALLOC(y+1,sizeof(P)); - mhp_rec(c,x,t,r2); - get_eg(&eg0); - for ( i = 0; i <= y; i++ ) { - mulp(CO,r[i],t,&s); addp(CO,s,r2[i],&r[i]); - } - get_eg(&eg1); add_eg(&eg_comp,&eg0,&eg1); -} - -void Pdp_monomial_hilbert_poincare(NODE arg,LIST *rp) -{ - LIST g,v; - VL vl; - int m,n,i; - VECT b,x; - NODE t,nd; - Z z; - P hp,tv,mt,t1,u,w; - DP *p; - P *plist,*r; - struct order_spec *spec; - struct oEGT eg0,eg1; - - if ( get_opt("hf",&val) && val ) hf = 1; - else hf = 0; - g = (LIST)ARG0(arg); v = (LIST)ARG1(arg); - pltovl(v,&vl); - m = length(BDY(g)); MKVECT(b,m); p = (DP *)BDY(b); - for ( t = BDY(g), i = 0; t; t = NEXT(t), i++ ) - ptod(CO,vl,(P)BDY(t),&p[i]); - n = length(BDY(v)); MKVECT(x,n); p = (DP *)BDY(x); - for ( t = BDY(v), i = 0; t; t = NEXT(t), i++ ) - ptod(CO,vl,(P)BDY(t),&p[i]); - create_order_spec(0,0,&spec); initd(spec); - /* create (1,1-t,...,(1-t)^n) */ - plist = (P *)MALLOC((n+1)*sizeof(P)); - /* t1 = 1-t */ - makevar("t",&tv); chsgnp(tv,&mt); addp(CO,mt,(P)ONE,&t1); - for ( plist[0] = (P)ONE, i = 1; i <= n; i++ ) - mulp(CO,plist[i-1],t1,&plist[i]); - r = (P *)CALLOC(n+1,sizeof(P)); - make_reduced(b); - mhp_rec(b,x,tv,r); - for ( hp = 0, i = 0; i <= n; i++ ) { - mulp(CO,plist[i],r[i],&u); addp(CO,u,hp,&w); hp = w; - } - UTOZ(n,z); - if ( !hf ) { - nd = mknode(2,hp,z); - MKLIST(*rp,nd); - } else { - P gcd,q; - int s; - Z qd; - ezgcdp(CO,hp,plist[n],&gcd); - if ( NUM(gcd) ) { - s = n; - q = hp; - } else { - s = n-ZTOS(DC(gcd)); - sdivp(CO,hp,plist[n-s],&q); - } - if ( NUM(q) ) qd = 0; - else qd = DEG(DC(q)); - nd = mknode(4,hp,z,q,qd); - MKLIST(*rp,nd); - } -} -#else - void dl_print(DL d,int n) { int i; @@ -599,7 +468,7 @@ void mhp_rec(VECT b,VECT x,P t,P *r) DL *p,*q; DL pi,xj,d1; VECT c; -struct oEGT eg0,eg1; + struct oEGT eg0,eg1; i_all++; n = b->len; nv = x->len; p = (DL *)BDY(b); @@ -652,6 +521,69 @@ struct oEGT eg0,eg1; // get_eg(&eg1); add_eg(&eg_comp,&eg0,&eg1); } +P mhp_rec_weight(VECT b,VECT x,P t,int *w) +{ + int n,i,j,k,l,i2,nv,len,td; + int *d; + Z wj; + P twj,tmp,tmp2,ret,qadd,qcolon; + DL *p,*q; + DL pi,xj,d1; + VECT c; + + i_all++; + n = b->len; nv = x->len; p = (DL *)BDY(b); + if ( !n ) { + // I=<0> => HP(t)=1/(1-t^w1)...(1-t^wn) => Q(t)=1 + return (P)ONE; + } + if ( n == 1 && p[0]->td == 0 ) { + // I=<1> => HP(t)=0 => Q(t)=0 + return 0; + } + for ( i = 0; i < n; i++ ) { + d = p[i]->d; + for ( td = 0, j = 0; j < nv; j++ ) td += d[j]; + if (td > 1 ) break; + } + if ( i == n ) { + // I= => Q(t)=(1-t^wi1)...(1-t^win) + for ( ret = (P)ONE, i = 0; i < n; i++ ) { + d = p[i]->d; + for ( j = 0; j < nv; j++ ) if ( d[j] ) break; + STOZ(w[j],wj); pwrp(CO,t,wj,&tmp); + subp(CO,(P)ONE,tmp,&tmp2); mulp(CO,ret,tmp2,&tmp); ret = tmp; + } + return ret; + } + for ( j = 0, d = p[i]->d; j < nv; j++ ) + if ( d[j] ) break; + xj = BDY(x)[j]; + MKVECT(c,n); q = (DL *)BDY(c); + for ( i = k = l = 0; i < n; i++ ) + if ( p[i]->d[j] ) { + pi = p[i]; + NEWDL(d1,nv); d1->td =pi->td - 1; + memcpy(d1->d,pi->d,nv*sizeof(int)); + d1->d[j]--; + p[k++] = d1; + } else + q[l++] = p[i]; + for ( i = k, i2 = 0; i2 < l; i++, i2++ ) + p[i] = q[i2]; + /* b=(b[0]/xj,...,b[k-1]/xj,b[k],...b[n-1]) where + b[0],...,b[k-1] are divisible by k */ + make_reduced2(b,k,nv); + qcolon = mhp_rec_weight(b,x,t,w); + /* c = (b[0],...,b[l-1],xj) */ + q[l] = xj; c->len = l+1; + qadd = mhp_rec_weight(c,x,t,w); + // Q(t)=Qadd+t^wj*Qcolon + STOZ(w[j],wj); pwrp(CO,t,wj,&twj); + mulp(CO,twj,qcolon,&tmp); addp(CO,qadd,tmp,&ret); + return ret; +} + /* (n+a)Cb as a polynomial of n; return (n+a)*...*(n+a-b+1) */ P binpoly(P n,int a,int b) @@ -668,17 +600,19 @@ P binpoly(P n,int a,int b) return r; } -void mhp_to_hf(VL vl,P hp,int n,P *plist,VECT *head,P *hf,Z *den) +void ibin(unsigned long int n,unsigned long int k,Z *r); + +void mhp_to_hf(VL vl,P hp,int n,P *plist,VECT *head,P *hf) { P tv,gcd,q,h,hphead,tt,ai,hpoly,nv,bp,w; - Z d; + Z d,z; DCP dc,topdc; VECT hfhead; int i,s,qd; if ( !hp ) { MKVECT(hfhead,0); *head = hfhead; - *hf = 0; *den = ONE; + *hf = 0; } else { makevar("t",&tv); ezgcdp(CO,hp,plist[n],&gcd); @@ -698,13 +632,12 @@ void mhp_to_hf(VL vl,P hp,int n,P *plist,VECT *head,P } *head = hfhead; *hf = 0; - *den = ONE; } else { if ( qd ) { topdc = 0; for ( i = 0; i < qd; i++ ) { NEWDC(dc); NEXT(dc) = topdc; - ibin(i+s-1,s-1,&COEF(dc)); + ibin(i+s-1,s-1,(Z *)&COEF(dc)); STOZ(i,d); DEG(dc) = d; topdc = dc; } @@ -724,8 +657,17 @@ void mhp_to_hf(VL vl,P hp,int n,P *plist,VECT *head,P addp(CO,hpoly,tt,&w); hpoly = w; } + if ( s > 2 ) { + factorialz(s-1,&z); + divsp(CO,hpoly,(P)z,&tt); hpoly = tt; + } *hf = hpoly; - factorialz(s-1,den); + for ( i = qd-1; i >= 0; i-- ) { + UTOZ(i,z); + substp(CO,hpoly,VR(nv),(P)z,&tt); + if ( cmpz((Z)tt,(Z)BDY(hfhead)[i]) ) break; + } + hfhead->len = i+1; } } } @@ -736,6 +678,7 @@ P *mhp_prep(int n,P *tv) { P *plist; P mt,t1; int i; + VECT list; plist = (P *)MALLOC((n+1)*sizeof(P)); /* t1 = 1-t */ @@ -756,20 +699,72 @@ P mhp_ctop(P *r,P *plist,int n) return hp; } +LIST dp_monomial_hilbert_poincare(VECT b,VECT x) +{ + int n; + P *r,*plist; + P tv; + P hp,hpoly; + VECT hfhead; + Z z; + NODE nd; + VECT vect; + LIST list; + + n = x->len; + plist = mhp_prep(n,&tv); + r = (P *)CALLOC(n+1,sizeof(P)); + make_reduced(b,n); + mhp_rec(b,x,tv,r); + hp = mhp_ctop(r,plist,n); + mhp_to_hf(CO,hp,n,plist,&hfhead,&hpoly); + UTOZ(n,z); + NEWVECT(vect); vect->len = n+1; BDY(vect) = (pointer)plist; + nd = mknode(5,hp,z,hfhead,hpoly,vect); + MKLIST(list,nd); + return list; +} + +LIST dp_monomial_hilbert_poincare_weight(VECT b,VECT x,int *w) +{ + int n,i; + NODE nd; + LIST list; + P tv,ret; + + n = x->len; + make_reduced(b,n); + makevar("t",&tv); + ret = mhp_rec_weight(b,x,tv,w); + nd = mknode(1,ret); + MKLIST(list,nd); + return list; +} + void Pdp_monomial_hilbert_poincare(NODE arg,LIST *rp) { LIST g,v; VL vl; - int m,n,i; - VECT b,x,hfhead; + int m,n,i,wlen; + VECT b,x,hfhead,prep; NODE t,nd; Z z,den; - P hp,tv,mt,t1,u,w,hpoly; + P hp,tv,mt,t1,u,hpoly; DP a; DL *p; - P *plist,*r; - Obj val; - + Obj val,ord,weight; + int *w; + struct order_spec *current_spec=0,*spec; + + weight = 0; + if ( current_option ) { + if ( peek_option(current_option,"ord",&ord) ) { + current_spec = dp_current_spec; + create_order_spec(0,ord,&spec); + initd(spec); + } + peek_option(current_option,"weight",&weight); + } i_simple = i_all = 0; g = (LIST)ARG0(arg); v = (LIST)ARG1(arg); pltovl(v,&vl); @@ -785,20 +780,92 @@ void Pdp_monomial_hilbert_poincare(NODE arg,LIST *rp) for ( t = BDY(v), i = 0; t; t = NEXT(t), i++ ) { ptod(CO,vl,(P)BDY(t),&a); p[i] = BDY(a)->dl; } + if ( weight ) { + wlen = length(BDY((LIST)weight)); + if ( n != wlen ) + error("dp_monomial_hilbert_poincare: inconsistent weight length"); + w = (int *)MALLOC(n*sizeof(int)); + for ( i = 0, nd = BDY((LIST)weight); i < n; i++, nd = NEXT(nd) ) + w[i] = ZTOS((Z)BDY(nd)); + } else if ( current_dl_weight_vector ) + w = current_dl_weight_vector; + else + w = 0; + if ( w ) { + *rp = dp_monomial_hilbert_poincare_weight(b,x,w); + } else { + *rp = dp_monomial_hilbert_poincare(b,x); + } + if ( current_spec ) + initd(current_spec); +} + +DL monomial_colon(DL a,DL b,int n) +{ + int i,d,td; + DL r; + + NEWDL(r,n); + td = 0; + for ( i = 0; i < n; i++ ) { + d = a->d[i]-b->d[i]; + r->d[i] = MAX(d,0); + td += r->d[i]; + } + r->td = td; + return r; +} + +// arguments : DPlist, Xlist, Mono, [HN(t),NV,Head,HP(n),Plist] +void Pdp_monomial_hilbert_poincare_incremental(NODE arg,LIST *rp) +{ + NODE g,data,data1,nd,t; + LIST list,list1; + DL new,dl; + int len,i,n; + Z dz; + DL *p; + P hn,hn1,newhn,tv,newhpoly,td,s; + VECT b,x,newhfhead; + P *plist; + Obj ord; + struct order_spec *current_spec=0,*spec; - r = (P *)CALLOC(n+1,sizeof(P)); + if ( current_option ) { + if ( peek_option(current_option,"ord",&ord) ) { + current_spec = dp_current_spec; + create_order_spec(0,ord,&spec); + initd(spec); + } + } + g = BDY((LIST)ARG0(arg)); new = BDY((DP)ARG1(arg))->dl; + data = BDY((LIST)ARG2(arg)); + hn = (P)ARG0(data); n = ZTOS((Z)ARG1(data)); + len = length(g); MKVECT(b,len); p = (DL *)BDY(b); + for ( t = g, i = 0; t; t = NEXT(t), i++ ) + p[i] = monomial_colon(BDY((DP)BDY(t))->dl,new,n); + MKVECT(x,n); + for ( i = 0; i < n; i++ ) { + NEWDL(dl,n); dl->d[i] = 1; dl->td = 1; BDY(x)[i] = dl; + } + // compute HP(I:new) + list1 = dp_monomial_hilbert_poincare(b,x); + data1 = BDY((LIST)list1); + hn1 = (P)ARG0(data1); + // HP(I+) = H(I)-t^d*H(I:new), d=tdeg(new) plist = mhp_prep(n,&tv); - make_reduced(b,n); - mhp_rec(b,x,tv,r); - hp = mhp_ctop(r,plist,n); - mhp_to_hf(CO,hp,n,plist,&hfhead,&hpoly,&den); - UTOZ(n,z); - nd = mknode(5,hp,z,hfhead,hpoly,den); - MKLIST(*rp,nd); + UTOZ(new->td,dz); + pwrp(CO,tv,dz,&td); + mulp(CO,hn1,td,&s); + subp(CO,hn,s,&newhn); + mhp_to_hf(CO,newhn,n,plist,&newhfhead,&newhpoly); + nd = mknode(5,newhn,ARG1(data),newhfhead,newhpoly,(VECT)ARG4(data)); + MKLIST(list,nd); + *rp = list; + if ( current_spec ) + initd(current_spec); } -#endif - void Pdp_compute_last_t(NODE arg,LIST *rp) { NODE g,gh,homo,n; @@ -955,6 +1022,19 @@ void Pdp_cont(NODE arg,Z *rp) dp_cont((DP)ARG0(arg),rp); } +void dpm_ptozp(DPM p,Z *cont,DPM *r); + +void Pdpm_remove_cont(NODE arg,LIST *rp) +{ + NODE nd; + Z cont; + DPM p; + + dpm_ptozp((DPM)ARG0(arg),&cont,&p); + nd = mknode(2,cont,p); + MKLIST(*rp,nd); +} + void Pdp_dtov(NODE arg,VECT *rp) { dp_dtov((DP)ARG0(arg),rp); @@ -1056,6 +1136,8 @@ void Pdp_nf_tab_f(NODE arg,DP *rp) dp_nf_tab_f((DP)ARG0(arg),(LIST *)BDY((VECT)ARG1(arg)),rp); } +extern int dpm_ordtype; + void Pdp_ord(NODE arg,Obj *rp) { struct order_spec *spec; @@ -1073,6 +1155,7 @@ void Pdp_ord(NODE arg,Obj *rp) else if ( !create_order_spec(0,(Obj)ARG0(arg),&spec) ) error("dp_ord : invalid order specification"); initd(spec); *rp = spec->obj; + if ( spec->id >= 256 ) dpm_ordtype = spec->module_ordtype; } } @@ -1223,7 +1306,7 @@ void Pdpm_ltod(NODE arg,DPM *rp) nd = BDY(f); len = length(nd); - for ( i = 0, t = nd, s = 0; i < len; i++, t = NEXT(t) ) { + for ( i = 1, t = nd, s = 0; i <= len; i++, t = NEXT(t) ) { ptod(CO,vl,(P)BDY(t),&d); dtodpm(d,i,&u); adddpm(CO,s,u,&w); s = w; @@ -1231,6 +1314,36 @@ void Pdpm_ltod(NODE arg,DPM *rp) *rp = s; } +// c*[monomial,i]+... -> c*<>+... + +void Pdpm_dptodpm(NODE arg,DPM *rp) +{ + DP p; + MP mp; + int pos,shift; + DMM m0,m; + + p = (DP)ARG0(arg); + pos = ZTOS((Z)ARG1(arg)); + if ( pos <= 0 ) + error("dpm_mtod : position must be positive"); + if ( !p ) *rp = 0; + else { + for ( m0 = 0, mp = BDY(p); mp; mp = NEXT(mp) ) { + NEXTDMM(m0,m); m->dl = mp->dl; m->c = mp->c; m->pos = pos; + } + if ( dp_current_spec->module_top_weight ) { + if ( pos > dp_current_spec->module_rank ) + error("dpm_dptodpm : inconsistent order spec"); + shift = dp_current_spec->module_top_weight[pos-1]; + m->dl->td += shift; + } else + shift = 0; + + MKDPM(p->nv,m0,*rp); (*rp)->sugar = p->sugar+shift; + } +} + void Pdpm_dtol(NODE arg,LIST *rp) { DPM a; @@ -1245,6 +1358,10 @@ void Pdpm_dtol(NODE arg,LIST *rp) Obj s; a = (DPM)ARG0(arg); + if ( !a ) { + MKLIST(*rp,0); + return; + } for ( vl = 0, nd = BDY((LIST)ARG1(arg)), nv = 0; nd; nd = NEXT(nd), nv++ ) { if ( !vl ) { NEWVL(vl); tvl = vl; @@ -1255,7 +1372,8 @@ void Pdpm_dtol(NODE arg,LIST *rp) } if ( vl ) NEXT(tvl) = 0; - n = ZTOS((Q)ARG2(arg)); + for ( t = BDY(a), n = 0; t; t = NEXT(t) ) + if ( t->pos > n ) n = t->pos; w = (MP *)CALLOC(n,sizeof(MP)); for ( t = BDY(a), len = 0; t; t = NEXT(t) ) len++; wa = (DMM *)MALLOC(len*sizeof(DMM)); @@ -1263,8 +1381,8 @@ void Pdpm_dtol(NODE arg,LIST *rp) for ( i = len-1; i >= 0; i-- ) { NEWMP(m); m->dl = wa[i]->dl; C(m) = C(wa[i]); pos = wa[i]->pos; - NEXT(m) = w[pos]; - w[pos] = m; + NEXT(m) = w[pos-1]; + w[pos-1] = m; } nd = 0; for ( i = n-1; i >= 0; i-- ) { @@ -1372,6 +1490,21 @@ void Pdp_mod(NODE arg,DP *rp) dp_mod(p,mod,subst,rp); } +void dpm_mod(DPM,int,DPM *); + +void Pdpm_mod(NODE arg,DPM *rp) +{ + DPM p; + int mod; + NODE subst; + + asir_assert(ARG0(arg),O_DP,"dp_mod"); + asir_assert(ARG1(arg),O_N,"dp_mod"); + p = (DPM)ARG0(arg); mod = ZTOS((Q)ARG1(arg)); + dpm_mod(p,mod,rp); +} + + void Pdp_rat(NODE arg,DP *rp) { asir_assert(ARG0(arg),O_DP,"dp_rat"); @@ -1380,11 +1513,74 @@ void Pdp_rat(NODE arg,DP *rp) extern int DP_Multiple; +int dp_iszp(DP); +int dpm_iszp(DPM); + +DP dptozdp(DP g) +{ + DP gz; + + if ( dp_iszp(g) ) + gz = g; + else + dp_ptozp(g,&gz); + return gz; +} + +VECT dpvtozdpv(VECT v) +{ + DP *ps,*psz; + int len,i; + VECT r; + + ps = (DP *)BDY(v); + len = v->len; + for ( i = 0; i < len; i++ ) + if ( !dp_iszp(ps[i]) ) break; + if ( i == len ) return v; + MKVECT(r,len); + psz = (DP *)BDY(r); + for ( i = 0; i < len; i++ ) + psz[i] = dptozdp(ps[i]); + return r; +} + +DPM dpmtozdpm(DPM g) +{ + DPM gz; + Z cont; + + if ( dpm_iszp(g) ) + gz = g; + else + dpm_ptozp(g,&cont,&gz); + return gz; +} + +VECT dpmvtozdpmv(VECT v) +{ + DPM *ps,*psz; + int len,i; + VECT r; + + ps = (DPM *)BDY(v); + len = v->len; + for ( i = 0; i < len; i++ ) + if ( !dpm_iszp(ps[i]) ) break; + if ( i == len ) return v; + MKVECT(r,len); + psz = (DPM *)BDY(r); + for ( i = 0; i < len; i++ ) + psz[i] = dpmtozdpm(ps[i]); + return r; +} + void Pdp_nf(NODE arg,DP *rp) { NODE b; DP *ps; DP g; + VECT zv; int full; do_weyl = 0; dp_fcoeffs = 0; @@ -1395,7 +1591,10 @@ void Pdp_nf(NODE arg,DP *rp) if ( !(g = (DP)ARG1(arg)) ) { *rp = 0; return; } - b = BDY((LIST)ARG0(arg)); ps = (DP *)BDY((VECT)ARG2(arg)); + b = BDY((LIST)ARG0(arg)); + zv = dpvtozdpv((VECT)ARG2(arg)); + g = dptozdp(g); + ps = (DP *)BDY(zv); full = (Q)ARG3(arg) ? 1 : 0; dp_nf_z(b,g,ps,full,DP_Multiple,rp); } @@ -1405,6 +1604,7 @@ void Pdp_weyl_nf(NODE arg,DP *rp) NODE b; DP *ps; DP g; + VECT zv; int full; asir_assert(ARG0(arg),O_LIST,"dp_weyl_nf"); @@ -1414,7 +1614,9 @@ void Pdp_weyl_nf(NODE arg,DP *rp) if ( !(g = (DP)ARG1(arg)) ) { *rp = 0; return; } - b = BDY((LIST)ARG0(arg)); ps = (DP *)BDY((VECT)ARG2(arg)); + b = BDY((LIST)ARG0(arg)); + zv = dpvtozdpv((VECT)ARG2(arg)); + g = dptozdp(g); full = (Q)ARG3(arg) ? 1 : 0; do_weyl = 1; dp_nf_z(b,g,ps,full,DP_Multiple,rp); @@ -1424,40 +1626,123 @@ void Pdp_weyl_nf(NODE arg,DP *rp) void Pdpm_nf(NODE arg,DPM *rp) { NODE b; - DPM *ps; + VECT ps; DPM g; - int full; + int ac,full; if ( !(g = (DPM)ARG1(arg)) ) { *rp = 0; return; } do_weyl = 0; dp_fcoeffs = 0; - asir_assert(ARG0(arg),O_LIST,"dpm_nf"); - asir_assert(ARG1(arg),O_DPM,"dpm_nf"); - asir_assert(ARG2(arg),O_VECT,"dpm_nf"); - asir_assert(ARG3(arg),O_N,"dpm_nf"); - b = BDY((LIST)ARG0(arg)); ps = (DPM *)BDY((VECT)ARG2(arg)); - full = (Q)ARG3(arg) ? 1 : 0; + ac = argc(arg); + if ( ac < 3 ) + error("dpm_nf: invalid arguments"); + else if ( ac == 3 ) { + asir_assert(ARG1(arg),O_VECT,"dpm_nf"); + b = 0; + g = dpmtozdpm((DPM)ARG0(arg)); + ps = dpmvtozdpmv((VECT)ARG1(arg)); + full = (Q)ARG2(arg) ? 1 : 0; + } else if ( ac == 4 ) { + asir_assert(ARG0(arg),O_LIST,"dpm_nf"); + asir_assert(ARG2(arg),O_VECT,"dpm_nf"); + b = BDY((LIST)ARG0(arg)); + g = dpmtozdpm((DPM)ARG1(arg)); + ps = dpmvtozdpmv((VECT)ARG2(arg)); + full = (Q)ARG3(arg) ? 1 : 0; + } dpm_nf_z(b,g,ps,full,DP_Multiple,rp); } +DP *dpm_nf_and_quotient(NODE b,DPM g,VECT ps,DPM *rp,P *dnp); +DPM dpm_nf_and_quotient2(NODE b,DPM g,VECT ps,DPM *rp,P *dnp); + +void Pdpm_nf_and_quotient(NODE arg,LIST *rp) +{ + NODE b; + VECT ps; + DPM g,nm; + P dn; + VECT quo; + NODE n; + int ac; + + do_weyl = 0; dp_fcoeffs = 0; + ac = argc(arg); + if ( ac < 2 ) + error("dpm_nf_and_quotient : invalid arguments"); + else if ( ac == 2 ) { + asir_assert(ARG1(arg),O_VECT,"dpm_nf_and_quotient"); + b = 0; g = (DPM)ARG0(arg); ps = (VECT)ARG1(arg); + } else if ( ac == 3 ) { + asir_assert(ARG0(arg),O_LIST,"dpm_nf_and_quotient"); + asir_assert(ARG2(arg),O_VECT,"dpm_nf_and_quotient"); + b = BDY((LIST)ARG0(arg)); g = (DPM)ARG1(arg); ps = (VECT)ARG2(arg); + } + NEWVECT(quo); quo->len = ps->len; + if ( g ) { + quo->body = (pointer *)dpm_nf_and_quotient(b,g,ps,&nm,&dn); + } else { + quo->body = (pointer *)MALLOC(quo->len*sizeof(pointer)); + nm = 0; dn = (P)ONE; + } + n = mknode(3,nm,dn,quo); + MKLIST(*rp,n); +} + +void Pdpm_nf_and_quotient2(NODE arg,LIST *rp) +{ + NODE b; + VECT ps; + DPM g,nm,q; + P dn; + NODE n; + int ac; + + do_weyl = 0; dp_fcoeffs = 0; + ac = argc(arg); + if ( ac < 2 ) + error("dpm_nf_and_quotient2 : invalid arguments"); + else if ( ac == 2 ) { + asir_assert(ARG1(arg),O_VECT,"dpm_nf_and_quotient2"); + b = 0; g = (DPM)ARG0(arg); ps = (VECT)ARG1(arg); + } else if ( ac == 3 ) { + asir_assert(ARG0(arg),O_LIST,"dpm_nf_and_quotient2"); + asir_assert(ARG2(arg),O_VECT,"dpm_nf_and_quotient2"); + b = BDY((LIST)ARG0(arg)); g = (DPM)ARG1(arg); ps = (VECT)ARG2(arg); + } + if ( g ) { + q = dpm_nf_and_quotient2(b,g,ps,&nm,&dn); + } else { + q = 0; nm = 0; dn = (P)ONE; + } + n = mknode(3,nm,dn,q); + MKLIST(*rp,n); +} + void Pdpm_weyl_nf(NODE arg,DPM *rp) { NODE b; - DPM *ps; + VECT ps; DPM g; - int full; + int ac,full; if ( !(g = (DPM)ARG1(arg)) ) { *rp = 0; return; } - asir_assert(ARG0(arg),O_LIST,"dpm_weyl_nf"); - asir_assert(ARG1(arg),O_DPM,"dpm_weyl_nf"); - asir_assert(ARG2(arg),O_VECT,"dpm_weyl_nf"); - asir_assert(ARG3(arg),O_N,"dpm_weyl_nf"); - b = BDY((LIST)ARG0(arg)); ps = (DPM *)BDY((VECT)ARG2(arg)); - full = (Q)ARG3(arg) ? 1 : 0; - do_weyl = 1; + do_weyl = 1; dp_fcoeffs = 0; + ac = argc(arg); + if ( ac < 3 ) + error("dpm_weyl_nf: invalid arguments"); + else if ( ac == 3 ) { + asir_assert(ARG1(arg),O_VECT,"dpm_nf"); + b = 0; g = (DPM)ARG0(arg); ps = (VECT)ARG1(arg); + } else if ( ac == 4 ) { + asir_assert(ARG0(arg),O_LIST,"dpm_weyl_nf"); + asir_assert(ARG2(arg),O_VECT,"dpm_weyl_nf"); + b = BDY((LIST)ARG0(arg)); g = (DPM)ARG1(arg); ps = (VECT)ARG2(arg); + full = (Q)ARG3(arg) ? 1 : 0; + } dpm_nf_z(b,g,ps,full,DP_Multiple,rp); do_weyl = 0; } @@ -1508,38 +1793,51 @@ void Pdp_weyl_nf_f(NODE arg,DP *rp) void Pdpm_nf_f(NODE arg,DPM *rp) { NODE b; - DPM *ps; + VECT ps; DPM g; - int full; + int ac,full; if ( !(g = (DPM)ARG1(arg)) ) { *rp = 0; return; } - asir_assert(ARG0(arg),O_LIST,"dpm_nf_f"); - asir_assert(ARG1(arg),O_DPM,"dpm_nf_f"); - asir_assert(ARG2(arg),O_VECT,"dpm_nf_f"); - asir_assert(ARG3(arg),O_N,"dpm_nf_f"); - b = BDY((LIST)ARG0(arg)); ps = (DPM *)BDY((VECT)ARG2(arg)); - full = (Q)ARG3(arg) ? 1 : 0; + ac = argc(arg); + if ( ac < 3 ) + error("dpm_nf_f: invalid arguments"); + else if ( ac == 3 ) { + asir_assert(ARG1(arg),O_VECT,"dpm_nf_f"); + b = 0; g = (DPM)ARG0(arg); ps = (VECT)ARG1(arg); + } else if ( ac == 4 ) { + asir_assert(ARG0(arg),O_LIST,"dpm_nf_f"); + asir_assert(ARG2(arg),O_VECT,"dpm_nf_f"); + b = BDY((LIST)ARG0(arg)); g = (DPM)ARG1(arg); ps = (VECT)ARG2(arg); + full = (Q)ARG3(arg) ? 1 : 0; + } + do_weyl = 0; dpm_nf_f(b,g,ps,full,rp); } void Pdpm_weyl_nf_f(NODE arg,DPM *rp) { NODE b; - DPM *ps; + VECT ps; DPM g; - int full; + int ac,full; if ( !(g = (DPM)ARG1(arg)) ) { *rp = 0; return; } - asir_assert(ARG0(arg),O_LIST,"dpm_weyl_nf_f"); - asir_assert(ARG1(arg),O_DP,"dpm_weyl_nf_f"); - asir_assert(ARG2(arg),O_VECT,"dpm_weyl_nf_f"); - asir_assert(ARG3(arg),O_N,"dpm_weyl_nf_f"); - b = BDY((LIST)ARG0(arg)); ps = (DPM *)BDY((VECT)ARG2(arg)); - full = (Q)ARG3(arg) ? 1 : 0; + ac = argc(arg); + if ( ac < 3 ) + error("dpm_weyl_nf_f: invalid arguments"); + else if ( ac == 3 ) { + asir_assert(ARG1(arg),O_VECT,"dpm_weyl_nf_f"); + b = 0; g = (DPM)ARG0(arg); ps = (VECT)ARG1(arg); + } else if ( ac == 4 ) { + asir_assert(ARG0(arg),O_LIST,"dpm_weyl_nf_f"); + asir_assert(ARG2(arg),O_VECT,"dpm_weyl_nf_f"); + b = BDY((LIST)ARG0(arg)); g = (DPM)ARG1(arg); ps = (VECT)ARG2(arg); + full = (Q)ARG3(arg) ? 1 : 0; + } do_weyl = 1; dpm_nf_f(b,g,ps,full,rp); do_weyl = 0; @@ -1579,9 +1877,11 @@ void Pdp_true_nf(NODE arg,LIST *rp) { NODE b,n; DP *ps; - DP g; - DP nm; - P dn; + DP g,gz; + DP nm,nm1; + P dn,dn1; + Z cont,cnm,cdn; + VECT zv; int full; do_weyl = 0; dp_fcoeffs = 0; @@ -1592,9 +1892,25 @@ void Pdp_true_nf(NODE arg,LIST *rp) if ( !(g = (DP)ARG1(arg)) ) { nm = 0; dn = (P)ONE; } else { - b = BDY((LIST)ARG0(arg)); ps = (DP *)BDY((VECT)ARG2(arg)); + b = BDY((LIST)ARG0(arg)); + zv = dpvtozdpv((VECT)ARG2(arg)); + ps = (DP *)BDY(zv); full = (Q)ARG3(arg) ? 1 : 0; - dp_true_nf(b,g,ps,full,&nm,&dn); + if ( dp_iszp(g) ) { + dp_true_nf(b,g,ps,full,&nm,&dn); + } else { + dp_ptozp3(g,&cont,&gz); + dp_true_nf(b,gz,ps,full,&nm1,&dn1); + if ( INT(cont) ) { + muldc(CO,nm1,(Obj)cont,&nm); + dn = dn1; + } else { + nmq((Q)cont,&cnm); + muldc(CO,nm1,(Obj)cnm,&nm); + dnq((Q)cont,&cdn); + mulp(CO,dn1,(P)cdn,&dn); + } + } } NEWNODE(n); BDY(n) = (pointer)nm; NEWNODE(NEXT(n)); BDY(NEXT(n)) = (pointer)dn; @@ -1922,6 +2238,93 @@ void Pdp_redble(NODE arg,Z *rp) *rp = 0; } +void Pdpm_redble(NODE arg,Z *rp) +{ + asir_assert(ARG0(arg),O_DPM,"dpm_redble"); + asir_assert(ARG1(arg),O_DPM,"dpm_redble"); + if ( dpm_redble((DPM)ARG0(arg),(DPM)ARG1(arg)) ) + *rp = ONE; + else + *rp = 0; +} + +void dpm_schreyer_base(LIST g,LIST *s); +void dpm_schreyer_base_zlist(LIST g,LIST *s); + +void Pdpm_schreyer_base(NODE arg,LIST *rp) +{ + asir_assert(ARG0(arg),O_LIST,"dpm_schreyer_base"); + dpm_schreyer_base_zlist((LIST)ARG0(arg),rp); +} + +void dpm_list_to_array(LIST g,VECT *psv,VECT *psiv); + +void Pdpm_list_to_array(NODE arg,LIST *rp) +{ + VECT psv,psiv; + NODE nd; + + asir_assert(ARG0(arg),O_LIST,"dpm_list_to_array"); + dpm_list_to_array((LIST)ARG0(arg),&psv,&psiv); + nd = mknode(2,psv,psiv); + MKLIST(*rp,nd); +} + +/* [quo,nf] = dpm_sp_nf(psv,psiv,i,j,top) */ +DPM dpm_sp_nf_zlist(VECT psv,VECT psiv,int i,int j,int top,DPM *nf); + +void Pdpm_sp_nf(NODE arg,LIST *rp) +{ + VECT psv,psiv; + DPM quo,nf; + Obj val; + int i,j,top; + NODE nd; + + asir_assert(ARG0(arg),O_VECT,"dpm_sp_nf"); psv = (VECT)ARG0(arg); + asir_assert(ARG1(arg),O_VECT,"dpm_sp_nf"); psiv = (VECT)ARG1(arg); + asir_assert(ARG2(arg),O_N,"dpm_sp_nf"); i = ZTOS((Q)ARG2(arg)); + asir_assert(ARG3(arg),O_N,"dpm_sp_nf"); j = ZTOS((Q)ARG3(arg)); + if ( get_opt("top",&val) && val ) + top = 1; + else + top = 0; + quo = dpm_sp_nf_zlist(psv,psiv,i,j,top,&nf); + nd = mknode(2,quo,nf); + MKLIST(*rp,nd); +} + +void dpm_insert_to_zlist(VECT psiv,int pos,int i); + +/* insert_to_zlist(indarray,dpm_hp(f),i) */ +void Pdpm_insert_to_zlist(NODE arg,VECT *rp) +{ + VECT psiv; + int i,pos; + + asir_assert(ARG0(arg),O_VECT,"dpm_insert_to_zlist"); psiv = (VECT)ARG0(arg); + asir_assert(ARG1(arg),O_N,"dpm_insert_to_zlist"); pos = ZTOS((Q)ARG1(arg)); + asir_assert(ARG2(arg),O_N,"dpm_insert_to_zlist"); i = ZTOS((Q)ARG2(arg)); + dpm_insert_to_zlist(psiv,pos,i); + *rp = psiv; +} + + +void dpm_simplify_syz(LIST m,LIST s,LIST *m1,LIST *s1,LIST *w1); + +void Pdpm_simplify_syz(NODE arg,LIST *rp) +{ + LIST s1,m1,w1; + NODE t; + + asir_assert(ARG0(arg),O_LIST,"dpm_simplify_syz"); + asir_assert(ARG1(arg),O_LIST,"dpm_simplify_syz"); + dpm_simplify_syz((LIST)ARG0(arg),(LIST)ARG1(arg),&s1,&m1,&w1); + t = mknode(3,s1,m1,w1); + MKLIST(*rp,t); +} + + void Pdp_red_mod(NODE arg,LIST *rp) { DP h,r; @@ -2084,25 +2487,47 @@ void Pdp_weyl_sp(NODE arg,DP *rp) do_weyl = 0; } -void Pdpm_sp(NODE arg,DPM *rp) +void Pdpm_sp(NODE arg,Obj *rp) { - DPM p1,p2; + DPM p1,p2,sp; + DP mul1,mul2; + Obj val; + NODE nd; + LIST l; do_weyl = 0; p1 = (DPM)ARG0(arg); p2 = (DPM)ARG1(arg); asir_assert(p1,O_DPM,"dpm_sp"); asir_assert(p2,O_DPM,"dpm_sp"); - dpm_sp(p1,p2,rp); + dpm_sp(p1,p2,&sp,&mul1,&mul2); + if ( get_opt("coef",&val) && val ) { + nd = mknode(3,sp,mul1,mul2); + MKLIST(l,nd); + *rp = (Obj)l; + } else { + *rp = (Obj)sp; + } } -void Pdpm_weyl_sp(NODE arg,DPM *rp) +void Pdpm_weyl_sp(NODE arg,Obj *rp) { - DPM p1,p2; + DPM p1,p2,sp; + DP mul1,mul2; + Obj val; + NODE nd; + LIST l; p1 = (DPM)ARG0(arg); p2 = (DPM)ARG1(arg); asir_assert(p1,O_DPM,"dpm_weyl_sp"); asir_assert(p2,O_DPM,"dpm_weyl_sp"); do_weyl = 1; - dpm_sp(p1,p2,rp); + dpm_sp(p1,p2,&sp,&mul1,&mul2); do_weyl = 0; + if ( get_opt("coef",&val) && val ) { + nd = mknode(3,sp,mul1,mul2); + MKLIST(l,nd); + *rp = (Obj)l; + } else { + *rp = (Obj)sp; + } } void Pdp_sp_mod(NODE arg,DP *rp) @@ -2183,6 +2608,17 @@ void Pdp_td(NODE arg,Z *rp) STOZ(BDY(p)->dl->td,*rp); } +void Pdpm_td(NODE arg,Z *rp) +{ + DPM p; + + p = (DPM)ARG0(arg); asir_assert(p,O_DPM,"dpm_td"); + if ( !p ) + *rp = 0; + else + STOZ(BDY(p)->dl->td,*rp); +} + void Pdp_sugar(NODE arg,Z *rp) { DP p; @@ -2423,6 +2859,22 @@ void Pdp_dehomo(NODE arg,DP *rp) dp_dehomo((DP)ARG0(arg),rp); } +void dpm_homo(DPM a,DPM *b); +void dpm_dehomo(DPM a,DPM *b); + +void Pdpm_homo(NODE arg,DPM *rp) +{ + asir_assert(ARG0(arg),O_DPM,"dpm_homo"); + dpm_homo((DPM)ARG0(arg),rp); +} + +void Pdpm_dehomo(NODE arg,DPM *rp) +{ + asir_assert(ARG0(arg),O_DPM,"dpm_dehomo"); + dpm_dehomo((DPM)ARG0(arg),rp); +} + + void Pdp_gr_flags(NODE arg,LIST *rp) { Obj name,value; @@ -2585,6 +3037,24 @@ void parse_gr_option(LIST f,NODE opt,LIST *v,Num *homo if ( !homo_is_set ) *homo = 0; } +int peek_option(NODE opt,char *find,Obj *retp) +{ + NODE t,p; + char *key; + Obj value; + + for ( t = opt; t; t = NEXT(t) ) { + p = BDY((LIST)BDY(t)); + key = BDY((STRING)BDY(p)); + value = (Obj)BDY(NEXT(p)); + if ( !strcmp(key,find) ) { + *retp = value; + return 1; + } + } + return 0; +} + void Pdp_gr_main(NODE arg,LIST *rp) { LIST f,v; @@ -2861,6 +3331,140 @@ void Pnd_gr(NODE arg,LIST *rp) nd_gr(f,v,m,homo,retdp,0,ord,rp); } +void nd_sba(LIST f,LIST v,int m,int homo,int retdp,int f4,struct order_spec *ord,LIST *rp); + +void Pnd_sba(NODE arg,LIST *rp) +{ + LIST f,v; + int m,homo,retdp,ac; + Obj val; + Z mq,z; + Num nhomo; + NODE node; + struct order_spec *ord,*current_spec; + + current_spec = dp_current_spec; + do_weyl = 0; + retdp = 0; + if ( (ac=argc(arg)) == 4 ) { + asir_assert(ARG0(arg),O_LIST,"nd_sba"); + asir_assert(ARG1(arg),O_LIST,"nd_sba"); + asir_assert(ARG2(arg),O_N,"nd_sba"); + f = (LIST)ARG0(arg); v = (LIST)ARG1(arg); + f = remove_zero_from_list(f); + if ( !BDY(f) ) { + *rp = f; return; + } + mq = (Z)ARG2(arg); + STOZ(0x40000000,z); + if ( cmpz(mq,z) >= 0 ) { + node = mknode(1,mq); + Psetmod_ff(node,&val); + m = -2; + } else + m = ZTOS(mq); + create_order_spec(0,ARG3(arg),&ord); + homo = 0; + if ( get_opt("homo",&val) && val ) homo = 1; + if ( get_opt("dp",&val) && val ) retdp = 1; + } else if ( ac == 1 ) { + f = (LIST)ARG0(arg); + parse_gr_option(f,current_option,&v,&nhomo,&m,&ord); + homo = ZTOS((Q)nhomo); + if ( get_opt("dp",&val) && val ) retdp = 1; + } else + error("nd_gr : invalid argument"); + nd_sba(f,v,m,homo,retdp,0,ord,rp); + initd(current_spec); +} + +void Pnd_weyl_sba(NODE arg,LIST *rp) +{ + LIST f,v; + int m,homo,retdp,ac; + Obj val; + Z mq,z; + Num nhomo; + NODE node; + struct order_spec *ord; + + do_weyl = 1; + retdp = 0; + if ( (ac=argc(arg)) == 4 ) { + asir_assert(ARG0(arg),O_LIST,"nd_sba"); + asir_assert(ARG1(arg),O_LIST,"nd_sba"); + asir_assert(ARG2(arg),O_N,"nd_sba"); + f = (LIST)ARG0(arg); v = (LIST)ARG1(arg); + f = remove_zero_from_list(f); + if ( !BDY(f) ) { + *rp = f; do_weyl = 0; return; + } + mq = (Z)ARG2(arg); + STOZ(0x40000000,z); + if ( cmpz(mq,z) >= 0 ) { + node = mknode(1,mq); + Psetmod_ff(node,&val); + m = -2; + } else + m = ZTOS(mq); + create_order_spec(0,ARG3(arg),&ord); + homo = 0; + if ( get_opt("homo",&val) && val ) homo = 1; + if ( get_opt("dp",&val) && val ) retdp = 1; + } else if ( ac == 1 ) { + f = (LIST)ARG0(arg); + parse_gr_option(f,current_option,&v,&nhomo,&m,&ord); + homo = ZTOS((Q)nhomo); + if ( get_opt("dp",&val) && val ) retdp = 1; + } else + error("nd_gr : invalid argument"); + nd_sba(f,v,m,homo,retdp,0,ord,rp); + do_weyl = 0; +} + +void Pnd_sba_f4(NODE arg,LIST *rp) +{ + LIST f,v; + int m,homo,retdp,ac; + Obj val; + Z mq,z; + Num nhomo; + NODE node; + struct order_spec *ord; + + do_weyl = 0; + retdp = 0; + if ( (ac=argc(arg)) == 4 ) { + asir_assert(ARG0(arg),O_LIST,"nd_sba"); + asir_assert(ARG1(arg),O_LIST,"nd_sba"); + asir_assert(ARG2(arg),O_N,"nd_sba"); + f = (LIST)ARG0(arg); v = (LIST)ARG1(arg); + f = remove_zero_from_list(f); + if ( !BDY(f) ) { + *rp = f; return; + } + mq = (Z)ARG2(arg); + STOZ(0x40000000,z); + if ( cmpz(mq,z) >= 0 ) { + node = mknode(1,mq); + Psetmod_ff(node,&val); + m = -2; + } else + m = ZTOS(mq); + create_order_spec(0,ARG3(arg),&ord); + homo = 0; + if ( get_opt("homo",&val) && val ) homo = 1; + if ( get_opt("dp",&val) && val ) retdp = 1; + } else if ( ac == 1 ) { + f = (LIST)ARG0(arg); + parse_gr_option(f,current_option,&v,&nhomo,&m,&ord); + homo = ZTOS((Q)nhomo); + if ( get_opt("dp",&val) && val ) retdp = 1; + } else + error("nd_gr : invalid argument"); + nd_sba(f,v,m,homo,retdp,1,ord,rp); +} + void Pnd_gr_postproc(NODE arg,LIST *rp) { LIST f,v; @@ -2972,6 +3576,8 @@ void Pnd_gr_trace(NODE arg,LIST *rp) { LIST f,v; int m,homo,ac; + Obj val; + int retdp; Num nhomo; struct order_spec *ord; @@ -2995,13 +3601,17 @@ void Pnd_gr_trace(NODE arg,LIST *rp) homo = ZTOS((Q)nhomo); } else error("nd_gr_trace : invalid argument"); - nd_gr_trace(f,v,m,homo,0,ord,rp); + retdp = 0; + if ( get_opt("dp",&val) && val ) retdp = 1; + nd_gr_trace(f,v,m,homo,retdp,0,ord,rp); } void Pnd_f4_trace(NODE arg,LIST *rp) { LIST f,v; int m,homo,ac; + int retdp; + Obj val; Num nhomo; struct order_spec *ord; @@ -3025,7 +3635,9 @@ void Pnd_f4_trace(NODE arg,LIST *rp) homo = ZTOS((Q)nhomo); } else error("nd_gr_trace : invalid argument"); - nd_gr_trace(f,v,m,homo,1,ord,rp); + retdp = 0; + if ( get_opt("dp",&val) && val ) retdp = 1; + nd_gr_trace(f,v,m,homo,retdp,1,ord,rp); } void Pnd_weyl_gr(NODE arg,LIST *rp) @@ -3066,7 +3678,8 @@ void Pnd_weyl_gr(NODE arg,LIST *rp) void Pnd_weyl_gr_trace(NODE arg,LIST *rp) { LIST f,v; - int m,homo,ac; + int m,homo,ac,retdp; + Obj val; Num nhomo; struct order_spec *ord; @@ -3090,7 +3703,9 @@ void Pnd_weyl_gr_trace(NODE arg,LIST *rp) homo = ZTOS((Q)nhomo); } else error("nd_weyl_gr_trace : invalid argument"); - nd_gr_trace(f,v,m,homo,0,ord,rp); + retdp = 0; + if ( get_opt("dp",&val) && val ) retdp = 1; + nd_gr_trace(f,v,m,homo,retdp,0,ord,rp); do_weyl = 0; } @@ -3909,25 +4524,89 @@ void Pdpv_ord(NODE arg,Obj *rp) *rp = dp_current_modspec->obj; } -extern int dpm_ispot; +extern int dpm_ordtype; +extern DMMstack dmm_stack; -void Pdpm_ord(NODE arg,LIST *rp) +void set_schreyer_order(LIST n); + +void Pdpm_set_schreyer(NODE arg,LIST *rp) { - Z q; - NODE nd; - struct order_spec *spec; - - if ( arg ) { - nd = BDY((LIST)ARG0(arg)); - if ( !create_order_spec(0,(Obj)ARG1(nd),&spec) ) - error("dpm_ord : invalid order specification"); - initdpm(spec,ZTOS((Q)ARG0(nd))); + if ( argc(arg) ) { + set_schreyer_order(ARG0(arg)?(LIST)ARG0(arg):0); } - STOZ(dpm_ispot,q); - nd = mknode(2,q,dp_current_spec->obj); - MKLIST(*rp,nd); + if ( dmm_stack ) + *rp = dmm_stack->obj; + else + *rp = 0; } +DMMstack_array Schreyer_Frame; +DMMstack_array dpm_schreyer_frame(NODE n,int lex); +void set_schreyer_level(DMMstack_array array,int level); + +void Pdpm_set_schreyer_level(NODE arg,Q *rp) +{ + set_schreyer_level(Schreyer_Frame,ZTOS((Q)ARG0(arg))); + *rp = (Q)ARG0(arg); +} + +void Pdpm_schreyer_frame(NODE arg,LIST *rp) +{ + DMMstack_array a; + DMMstack *body; + DMM *in,*sum; + DPM f,s; + NODE b,b1,nd; + LIST l; + VECT v; + Z lev,deg,ind; + int len,i,nv,rank,j,lex; + NODE tt,p; + char *key; + Obj value; + + lex = 0; + if ( current_option ) { + for ( tt = current_option; tt; tt = NEXT(tt) ) { + p = BDY((LIST)BDY(tt)); + key = BDY((STRING)BDY(p)); + value = (Obj)BDY(NEXT(p)); + if ( !strcmp(key,"lex") ) + lex = value!=0?1:0; + else { + error("dpm_schreyer_frame: unknown option."); + } + } + } + Schreyer_Frame = a = dpm_schreyer_frame(BDY((LIST)ARG0(arg)),lex); + len = a->len; + body = a->body; + /* XXX */ + nv = ((DPM)BDY(BDY((LIST)body[0]->obj)))->nv; + b = 0; + for ( i = 0; i < len; i++ ) { + rank = body[i]->rank; + in = body[i]->in; + sum = body[i]->sum; + MKVECT(v,rank+1); + STOZ(i+1,lev); + for ( j = 1; j <= rank; j++ ) { + MKDPM(nv,in[j],f); f->sugar = in[j]->dl->td; + MKDPM(nv,sum[j],s);s->sugar = sum[j]->dl->td; + STOZ(s->sugar,deg); + STOZ(j,ind); + nd = mknode(5,f,s,ind,lev,deg); + MKLIST(l,nd); + BDY(v)[j] = (pointer)l; + } + MKNODE(b1,(pointer)v,b); + b = b1; + } + MKLIST(l,b); + *rp = l; +} + + void Pdpm_hm(NODE arg,DPM *rp) { DPM p; @@ -3940,20 +4619,102 @@ void Pdpm_ht(NODE arg,DPM *rp) { DPM p; - p = (DPM)ARG0(arg); asir_assert(p,O_DPM,"dp_ht"); + p = (DPM)ARG0(arg); asir_assert(p,O_DPM,"dpm_ht"); dpm_ht(p,rp); } -void Pdpm_hc(NODE arg,Obj *rp) +void dpm_rest(DPM p,DPM *r); + +void Pdpm_rest(NODE arg,DPM *rp) { + DPM p; + + p = (DPM)ARG0(arg); asir_assert(p,O_DPM,"dpm_ht"); + dpm_rest(p,rp); +} + + +void Pdpm_hp(NODE arg,Z *rp) +{ + DPM p; + int pos; + + p = (DPM)ARG0(arg); asir_assert(p,O_DPM,"dpm_ht"); + pos = BDY(p)->pos; + STOZ(pos,*rp); +} + +void dpm_shift(DPM p,int s,DPM *rp); + +void Pdpm_shift(NODE arg,DPM *rp) +{ + DPM p; + int s; + + p = (DPM)ARG0(arg); asir_assert(p,O_DPM,"dpm_shift"); + s = ZTOS((Z)ARG1(arg)); + dpm_shift(p,s,rp); +} + +void dpm_sort(DPM p,DPM *rp); + +void Pdpm_sort(NODE arg,DPM *rp) +{ + DPM p; + int s; + + p = (DPM)ARG0(arg); + if ( !p ) *rp = 0; + else dpm_sort(p,rp); +} + +void dpm_split(DPM p,int s,DPM *up,DPM *lo); +void dpm_extract(DPM p,int s,DP *r); + +void Pdpm_split(NODE arg,LIST *rp) +{ + DPM p,up,lo; + int s; + NODE nd; + + p = (DPM)ARG0(arg); + s = ZTOS((Z)ARG1(arg)); + dpm_split(p,s,&up,&lo); + nd = mknode(2,up,lo); + MKLIST(*rp,nd); +} + +void Pdpm_extract(NODE arg,DP *rp) +{ + DPM p; + int s; + + p = (DPM)ARG0(arg); + s = ZTOS((Z)ARG1(arg)); + dpm_extract(p,s,rp); +} + + +void Pdpm_hc(NODE arg,DP *rp) +{ + DPM p; + DP d; + MP m; + asir_assert(ARG0(arg),O_DPM,"dpm_hc"); if ( !ARG0(arg) ) *rp = 0; - else - *rp = BDY((DPM)ARG0(arg))->c; + else { + p = (DPM)ARG0(arg); + NEWMP(m); + m->dl = BDY(p)->dl; + m->c = BDY(p)->c; + NEXT(m) = 0; + MKDP(NV(p),m,d); d->sugar = p->sugar; + *rp = d; + } } - void Pdpv_ht(NODE arg,LIST *rp) { NODE n; @@ -4053,6 +4814,10 @@ int dpv_hp(DPV p) case ORD_LEX: for ( i = 0; i < len; i++ ) if ( e[i] ) return i; + return -1; + break; + default: + error("dpv_hp : unsupported term ordering"); return -1; break; }