===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2018/builtin/dp.c,v
retrieving revision 1.29
retrieving revision 1.32
diff -u -p -r1.29 -r1.32
--- OpenXM_contrib2/asir2018/builtin/dp.c	2021/02/28 02:33:15	1.29
+++ 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.28 2021/02/18 05:35:01 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"
@@ -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);
@@ -520,6 +521,69 @@ void mhp_rec(VECT b,VECT x,P t,P *r)
 // 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)
@@ -635,10 +699,10 @@ P mhp_ctop(P *r,P *plist,int n)
   return hp;
 }
 
-LIST dp_monomial_hilbert_poincare(VECT b,VECT x,P *plist)
+LIST dp_monomial_hilbert_poincare(VECT b,VECT x)
 {
   int n;
-  P *r;
+  P *r,*plist;
   P tv;
   P hp,hpoly;
   VECT hfhead;
@@ -648,9 +712,9 @@ LIST dp_monomial_hilbert_poincare(VECT b,VECT x,P *pli
   LIST list;
 
   n = x->len;
+  plist = mhp_prep(n,&tv);
   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);
@@ -661,20 +725,46 @@ LIST dp_monomial_hilbert_poincare(VECT b,VECT x,P *pli
   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;
+  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;
-  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);
@@ -690,8 +780,24 @@ 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;
   }
-  plist = mhp_prep(n,&tv);
-  *rp = dp_monomial_hilbert_poincare(b,x,plist);
+  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)
@@ -722,11 +828,19 @@ void Pdp_monomial_hilbert_poincare_incremental(NODE ar
   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);
@@ -735,11 +849,12 @@ void Pdp_monomial_hilbert_poincare_incremental(NODE ar
     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);
+  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)
-  makevar("t",&tv); UTOZ(new->td,dz);
+  plist = mhp_prep(n,&tv);
+  UTOZ(new->td,dz);
   pwrp(CO,tv,dz,&td);
   mulp(CO,hn1,td,&s);
   subp(CO,hn,s,&newhn);
@@ -747,6 +862,8 @@ void Pdp_monomial_hilbert_poincare_incremental(NODE ar
   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)
@@ -1396,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;
@@ -1411,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);
 }
@@ -1421,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");
@@ -1430,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);
@@ -1453,11 +1639,16 @@ void Pdpm_nf(NODE arg,DPM *rp)
     error("dpm_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);
+    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 = (DPM)ARG1(arg); ps = (VECT)ARG2(arg);
+    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);
@@ -1686,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;
@@ -1699,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;
@@ -2826,6 +3035,24 @@ void parse_gr_option(LIST f,NODE opt,LIST *v,Num *homo
   if ( !ord_is_set ) create_order_spec(0,0,ord);
   if ( !modular_is_set ) *modular = 0;
   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)