===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2018/builtin/dp.c,v
retrieving revision 1.25
retrieving revision 1.29
diff -u -p -r1.25 -r1.29
--- OpenXM_contrib2/asir2018/builtin/dp.c	2020/06/30 01:52:17	1.25
+++ OpenXM_contrib2/asir2018/builtin/dp.c	2021/02/28 02:33:15	1.29
@@ -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.24 2020/06/19 10:18:13 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2018/builtin/dp.c,v 1.28 2021/02/18 05:35:01 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();
@@ -114,6 +114,7 @@ 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();
@@ -196,6 +197,7 @@ struct ftab dp_tab[] = {
   {"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},
@@ -247,6 +249,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},
@@ -291,6 +294,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},
@@ -463,7 +467,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);
@@ -610,6 +614,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 */
@@ -630,18 +635,44 @@ 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;
+  P *plist;
   Obj val;
 
   i_simple = i_all = 0;
@@ -659,18 +690,65 @@ 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);
 }
 
+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;
+
+  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;
+}
+
 void Pdp_compute_last_t(NODE arg,LIST *rp)
 {
   NODE g,gh,homo,n;
@@ -1295,6 +1373,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");
@@ -3021,8 +3114,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 ) {
@@ -3054,6 +3148,51 @@ 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)
+{
+  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)