===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2018/builtin/dp.c,v
retrieving revision 1.23
retrieving revision 1.31
diff -u -p -r1.23 -r1.31
--- OpenXM_contrib2/asir2018/builtin/dp.c	2020/02/11 01:43:57	1.23
+++ OpenXM_contrib2/asir2018/builtin/dp.c	2021/12/05 22:41:03	1.31
@@ -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.22 2020/01/09 01:47:40 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2018/builtin/dp.c,v 1.30 2021/03/10 06:36:20 noro Exp $
 */
 #include "ca.h"
 #include "base.h"
@@ -61,7 +61,7 @@ extern int nd_rref2;
 
 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,7 +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();
+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();
@@ -113,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();
@@ -131,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);
@@ -194,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},
@@ -244,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},
@@ -288,6 +295,7 @@ struct ftab dp_supp_tab[] = {
   {"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},
@@ -460,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);
@@ -513,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=<xi1,...,xin> => 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)
@@ -607,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 */
@@ -627,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);
@@ -656,16 +780,90 @@ 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+<new>) = 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);
-  UTOZ(n,z);
-  nd = mknode(4,hp,z,hfhead,hpoly);
-  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);
 }
 
 void Pdp_compute_last_t(NODE arg,LIST *rp)
@@ -1292,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");
@@ -2732,6 +2945,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;
@@ -3006,6 +3237,140 @@ void Pnd_gr(NODE arg,LIST *rp)
   } else
     error("nd_gr : invalid argument");
   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)