===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2018/builtin/dp.c,v
retrieving revision 1.27
retrieving revision 1.30
diff -u -p -r1.27 -r1.30
--- OpenXM_contrib2/asir2018/builtin/dp.c	2021/01/25 00:39:51	1.27
+++ OpenXM_contrib2/asir2018/builtin/dp.c	2021/03/10 06:36:20	1.30
@@ -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.26 2020/10/26 02:41:05 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2018/builtin/dp.c,v 1.29 2021/02/28 02:33:15 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();
@@ -133,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);
@@ -249,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},
@@ -466,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);
@@ -613,6 +615,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 */
@@ -633,20 +636,54 @@ P mhp_ctop(P *r,P *plist,int n)
   return hp;
 }
 
+LIST dp_monomial_hilbert_poincare(VECT b,VECT x,P *plist)
+{
+  int n;
+  P *r;
+  P tv;
+  P hp,hpoly;
+  VECT hfhead;
+  Z z;
+  NODE nd;
+  VECT vect;
+  LIST list;
+
+  n = x->len;
+  r = (P *)CALLOC(n+1,sizeof(P));
+  make_reduced(b,n);
+  makevar("t",&tv);
+  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;
+}
+
 void Pdp_monomial_hilbert_poincare(NODE arg,LIST *rp)
 {
   LIST g,v;
   VL vl;
   int m,n,i;
-  VECT b,x,hfhead;
+  VECT b,x,hfhead,prep;
   NODE t,nd;
   Z z,den;
   P hp,tv,mt,t1,u,w,hpoly;
   DP a;
   DL *p;
-  P *plist,*r;
-  Obj val;
-
+  P *plist;
+  Obj val,ord;
+  struct order_spec *current_spec=0,*spec;
+  
+  if ( current_option ) {
+    if ( peek_option(current_option,"ord",&ord) ) {
+      current_spec = dp_current_spec;
+      create_order_spec(0,ord,&spec);
+      initd(spec);
+    }
+  }
   i_simple = i_all = 0;
   g = (LIST)ARG0(arg); v = (LIST)ARG1(arg);
   pltovl(v,&vl);
@@ -662,18 +699,78 @@ 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;
   }
-  
-  r = (P *)CALLOC(n+1,sizeof(P));
   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);
+  *rp = dp_monomial_hilbert_poincare(b,x,plist);
+  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;
+  
+  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)); 
+  plist = (P *)BDY((VECT)ARG4(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,plist);
+  data1 = BDY((LIST)list1);
+  hn1 = (P)ARG0(data1);
+  // HP(I+<new>) = H(I)-t^d*H(I:new), d=tdeg(new)
+  makevar("t",&tv); 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)
 {
   NODE g,gh,homo,n;
@@ -2753,6 +2850,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;
@@ -3039,8 +3154,9 @@ void Pnd_sba(NODE arg,LIST *rp)
   Z mq,z;
   Num nhomo;
   NODE node;
-  struct order_spec *ord;
+  struct order_spec *ord,*current_spec;
 
+  current_spec = dp_current_spec;
   do_weyl = 0;
   retdp = 0;
   if ( (ac=argc(arg)) == 4 ) {
@@ -3072,6 +3188,7 @@ void Pnd_sba(NODE arg,LIST *rp)
   } 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)