===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2018/builtin/dp.c,v
retrieving revision 1.22
retrieving revision 1.30
diff -u -p -r1.22 -r1.30
--- OpenXM_contrib2/asir2018/builtin/dp.c	2020/01/09 01:47:40	1.22
+++ 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.21 2019/12/27 08:13:59 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();
@@ -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);
@@ -607,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 */
@@ -627,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);
@@ -656,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;
@@ -1292,6 +1395,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 +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;
@@ -3008,6 +3144,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;
@@ -4084,7 +4354,7 @@ void Pdpm_set_schreyer(NODE arg,LIST *rp)
 }
 
 DMMstack_array Schreyer_Frame;
-DMMstack_array dpm_schreyer_frame(NODE n);
+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)
@@ -4093,10 +4363,6 @@ void Pdpm_set_schreyer_level(NODE arg,Q *rp)
   *rp = (Q)ARG0(arg);
 }
 
-DPM dmmtodpm(DMM d)
-{
-}
-
 void Pdpm_schreyer_frame(NODE arg,LIST *rp)
 {
   DMMstack_array a;
@@ -4107,9 +4373,25 @@ void Pdpm_schreyer_frame(NODE arg,LIST *rp)
   LIST l;
   VECT v;
   Z lev,deg,ind;
-  int len,i,nv,rank,j;
+  int len,i,nv,rank,j,lex;
+  NODE tt,p;
+  char *key;
+  Obj value;
 
-  Schreyer_Frame = a = dpm_schreyer_frame(BDY((LIST)ARG0(arg)));
+  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 */