===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2018/engine/dist.c,v
retrieving revision 1.3
retrieving revision 1.9
diff -u -p -r1.3 -r1.9
--- OpenXM_contrib2/asir2018/engine/dist.c	2019/08/21 00:37:47	1.3
+++ OpenXM_contrib2/asir2018/engine/dist.c	2019/09/19 06:29:47	1.9
@@ -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/engine/dist.c,v 1.2 2018/09/28 08:20:28 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2018/engine/dist.c,v 1.8 2019/09/13 02:04:42 noro Exp $
 */
 #include "ca.h"
 
@@ -333,7 +333,7 @@ void nodetod(NODE node,DP *dp)
 void nodetodpm(NODE node,Obj pos,DPM *dp)
 {
   NODE t;
-  int len,i,td;
+  int len,i,td,p;
   Q e;
   DL d;
   DMM m;
@@ -352,7 +352,13 @@ void nodetodpm(NODE node,Obj pos,DPM *dp)
     }
   }
   d->td = td;
-  NEWDMM(m); m->dl = d; m->pos = ZTOS((Q)pos); C(m) = (Obj)ONE; NEXT(m) = 0;
+  p = ZTOS((Q)pos);
+  if ( dp_current_spec->module_rank ) {
+    if ( p > dp_current_spec->module_rank )
+      error("nodetodpm : inconsistent order spec");
+    d->td += dp_current_spec->module_top_weight[p-1];
+  }
+  NEWDMM(m); m->dl = d; m->pos = p; C(m) = (Obj)ONE; NEXT(m) = 0;
   MKDPM(len,m,u); u->sugar = td; *dp = u;
 }
 
@@ -360,12 +366,20 @@ void dtodpm(DP d,int pos,DPM *dp)
 {
   DMM mr0,mr;
   MP m;
+  int shift;
 
   if ( !d ) *dp = 0;
   else {
+    shift = 0;
+    if ( dp_current_spec->module_rank ) {
+      if ( pos > dp_current_spec->module_rank )
+        error("nodetodpm : inconsistent order spec");
+      shift = dp_current_spec->module_top_weight[pos-1];
+    } 
     for ( m = BDY(d), mr0 = 0; m; m = NEXT(m) ) {
       NEXTDMM(mr0,mr);
       mr->dl = m->dl;
+      mr->dl->td += shift;
       mr->pos = pos;
       C(mr) = C(m);
     }
@@ -2635,73 +2649,218 @@ NBP harmonic_mul_nbm(NBM a,NBM b)
 /* DPM functions */
 
 DMMstack dmm_stack;
+int dpm_base_ordtype;;
 
-// data=[Ink,...,In0,Base]
-// Ini = a list of module monomials
-// Base=an order spec for polynomial ring or module
-void set_schreyer_order(NODE data)
+DMMstack push_schreyer_order(LIST data,DMMstack stack)
 {
   DMMstack t;
   int len,i;
+  NODE in,t1;
+  LIST l;
+
+  /* data = [DPM,...,DPM] */
+  if ( !dmm_stack && ( !dp_current_spec || dp_current_spec->id < 256 ) )
+      error("push_schreyer_order : base module order is not set");
+  in = BDY(data);
+  len = length(in);
+  NEWDMMstack(t);
+  t->rank = len;
+  t->in = (DMM *)MALLOC((len+1)*sizeof(DMM));
+  if ( stack ) {
+    MKNODE(t1,data,BDY(stack->obj)); MKLIST(l,t1); t->obj = l;
+  } else {
+    MKNODE(t1,data,0); MKLIST(l,t1); t->obj = l;
+  }
+  for ( i = 1; i <= len; i++, in = NEXT(in) ) {
+    t->in[i] = BDY((DPM)BDY(in));
+  }
+  t->next = stack;;
+  dpm_ordtype = 3;
+  return t;
+}
+
+// data=[Ink,...,In0]
+// Ini = a list of module monomials
+
+void set_schreyer_order(LIST data)
+{
   NODE in;
-  struct order_spec *base;
+  LIST *w;
+  int i,len;
 
   if ( !data ) {
     dmm_stack = 0;
     if ( dp_current_spec && dp_current_spec->id >= 256 )
-      dpm_ordtype = dp_current_spec->ispot;
+      dpm_ordtype = dp_current_spec->module_ordtype;
     else 
       dpm_ordtype = 0;
     return;
-  } else if ( NEXT(data) == 0 ) {
-    create_order_spec(0,BDY(data),&base);
-    NEWDMMstack(t);
-    t->in = 0;
-    t->rank = 0;
-    t->ordtype = base->ispot;
-    t->next = 0;
-    dmm_stack = t;
-    dpm_ordtype = 2;
   } else {
-    set_schreyer_order(NEXT(data));
-    in = BDY((LIST)BDY(data));
+    if ( !dp_current_spec || dp_current_spec->id < 256 )
+      error("set_schreyer_order : base module order is not set");
+    dmm_stack = 0;
+    dpm_base_ordtype = dp_current_spec->module_ordtype;
+    in = BDY(data);
     len = length(in);
-    NEWDMMstack(t);
-    t->in = 0;
-    t->rank = len;
-    t->in = (DMM *)MALLOC((len+1)*sizeof(DMM));
-    t->ordtype = 0;
-    t->next = dmm_stack;
-    dmm_stack = t; 
-    for ( i = 1; i <= len; i++, in = NEXT(in) ) {
-      t->in[i] = BDY((DPM)BDY(in));
+    w = (LIST *)MALLOC(len*sizeof(LIST));
+    for ( i = 0; i < len; i++, in = NEXT(in) ) w[i] = (LIST)BDY(in);
+    for ( i = len-1; i >= 0; i-- ) dmm_stack = push_schreyer_order(w[i],dmm_stack);
+  }
+}
+
+// construct a base of syz(g) 
+// assuming the schrerer order is properly set
+
+DP dpm_sp_hm(DPM p1,DPM p2);
+void dpm_sp(DPM p1,DPM p2,DPM *sp,DP *t1,DP *t2);
+DP *dpm_nf_and_quotient(NODE b,DPM sp,VECT psv,DPM *nf,P *dn);
+void dpm_sort(DPM p,DPM *r);
+
+extern int DP_Multiple;
+
+void dpm_nf_z(NODE b,DPM g,VECT psv,int full,int multiple,DPM *rp);
+NODE dpm_sort_list(NODE l);
+void dpm_ptozp(DPM p,Z *cont,DPM *r);
+
+NODE dpm_reduceall(NODE in)
+{
+  int n,i;
+  VECT psv;
+  DPM *ps;
+  NODE t,t1;
+  DPM g,r;
+  Z cont;
+
+  n = length(in);
+  MKVECT(psv,n);
+  ps = (DPM *)BDY(psv);
+  for ( i = 0, t = in; i < n; i++, t = NEXT(t) ) ps[i] = BDY(t);
+  for ( i = 0; i < n; i++ ) {
+    g = ps[i]; ps[i] = 0;
+//    dpm_nf_z(0,g,psv,1,DP_Multiple,&ps[i]);
+    dpm_nf_z(0,g,psv,1,0,&ps[i]);
+  }
+  t = 0;
+  for ( i = n-1; i >= 0; i-- ) {
+    dpm_ptozp(ps[i],&cont,&r);
+    MKNODE(t1,r,t); t = t1;
+  }
+  return t;
+}
+
+struct oEGT egra;
+
+void dpm_ht(DPM d,DPM *r);
+
+void dpm_schreyer_base(LIST g,LIST *s)
+{
+  NODE nd,t0,t,b0,b;
+  int n,i,j,k,nv;
+  LIST l;
+  Z cont;
+  P dn,c;
+  DP h,t1,t2;
+  MP d;
+  DMM r0,r;
+  DPM sp,nf,dpm;
+  DPM *ps;
+  VECT psv;
+  DP **m,*quo;
+  struct oEGT eg0,eg1;
+  extern struct oEGT egred;
+
+  init_eg(&egra);
+  nd = BDY(g);
+  n = length(nd);
+  MKVECT(psv,n);
+  ps = (DPM *)BDY(psv);
+  for ( i = 0, t = nd; i < n; i++, t = NEXT(t) ) ps[i] = (DPM)BDY(t);
+  nv = ps[0]->nv;
+  m = (DP **)almat_pointer(n,n);
+  b0 = 0;
+  for ( i = 0; i < n; i++ ) {
+    // sp(ps[i],ps[j]) = ti*ps[i]-tj*ps[j] => m[i][j] = ti
+    for ( j = i+1; j < n; j++ ) m[i][j] = dpm_sp_hm(ps[i],ps[j]);
+    for ( j = i+1; j < n; j++ ) {
+      if ( !m[i][j] ) continue;
+      for ( h = m[i][j], k = i+1; k < n; k++ )
+        if ( k != j && m[i][k] && dp_redble(m[i][k],h) ) m[i][k] = 0;
     }
+    for ( j = i+1; j < n; j++ ) {
+      if ( m[i][j] ) {
+        dpm_sp(ps[i],ps[j],&sp,&t1,&t2);
+        quo = dpm_nf_and_quotient(0,sp,psv,&nf,&dn);
+        if ( nf ) 
+          error("dpm_schreyer_base : cannot happen");
+        NEWDMM(r0); r = r0; 
+        mulp(CO,(P)BDY(t1)->c,dn,(P *)&r->c); r->pos = i+1; r->dl = BDY(t1)->dl;
+        NEWDMM(NEXT(r)); r=NEXT(r);
+        mulp(CO,(P)BDY(t2)->c,dn,&c); chsgnp(c,(P *)&r->c); r->pos = j+1; r->dl = BDY(t2)->dl;
+        if ( quo ) {
+          for ( k = 0; k < n; k++ ) {
+            if ( !quo[k] ) continue;
+            for ( d = BDY(quo[k]); d; d = NEXT(d) ) {
+              NEXTDMM(r0,r); chsgnp((P)d->c,(P *)&r->c); r->pos = k+1; r->dl = d->dl;
+            }
+          }
+        }
+        NEXT(r) = 0;
+        MKDPM(nv,r0,dpm); // XXX : sugar is not set
+        NEXTNODE(b0,b);
+        BDY(b) = (pointer)dpm;
+      }
+    }
+    if ( b0 ) NEXT(b) = 0;
   }
+  for ( t0 = t, nd = BDY(g); nd; nd = NEXT(nd) ) {
+    dpm_ht((DPM)BDY(nd),&dpm); NEXTNODE(t0,t); BDY(t) = (pointer)dpm;
+  }
+  if ( t0 ) NEXT(t) = 0;
+  MKLIST(l,t0);
+  dmm_stack = push_schreyer_order(l,dmm_stack);
+  for ( t = b0; t; t = NEXT(t) ) {
+    dpm_sort((DPM)BDY(t),&dpm); 
+    BDY(t) = (pointer)dpm;
+  }
+  b0 = dpm_sort_list(b0);
+// get_eg(&eg0);
+  b0 = dpm_reduceall(b0);
+// get_eg(&eg1); add_eg(&egra,&eg0,&eg1); print_eg("RA",&egra);
+  MKLIST(*s,b0);
+//  print_eg("red",&egred); printf("\n");
 }
 
 int compdmm_schreyer(int n,DMM m1,DMM m2)
 {
-   DL d1,d2;
    int pos1,pos2,t;
    DMM *in;
    DMMstack s;
+   static DL d1=0,d2=0;
+   static int dlen=0;
 
-   NEWDL(d1,n); _copydl(n,m1->dl,d1); pos1 = m1->pos;
-   NEWDL(d2,n); _copydl(n,m2->dl,d2); pos2 = m2->pos;
-   for ( s = dmm_stack; s->in; s = NEXT(s) ) {
+   pos1 = m1->pos; pos2 = m2->pos;
+   if ( pos1 == pos2 ) return (*cmpdl)(n,m1->dl,m2->dl);
+   if ( n > dlen ) {
+     NEWDL(d1,n); NEWDL(d2,n); dlen = n;
+   }
+   _copydl(n,m1->dl,d1); 
+   _copydl(n,m2->dl,d2); 
+   for ( s = dmm_stack; s; s = NEXT(s) ) {
      in = s->in;
      _addtodl(n,in[pos1]->dl,d1);
      _addtodl(n,in[pos2]->dl,d2);
-     if ( _eqdl(n,d1,d2) && in[pos1]->pos == in[pos2]->pos ) {
+     if ( in[pos1]->pos == in[pos2]->pos && _eqdl(n,d1,d2)) {
        if ( pos1 < pos2 ) return 1;
        else if ( pos1 > pos2 ) return -1;
        else return 0;
      }
      pos1 = in[pos1]->pos;
      pos2 = in[pos2]->pos;
+     if ( pos1 == pos2 ) return (*cmpdl)(n,d1,d2);
    }
    // comparison by the bottom order
-  if ( s->ordtype == 1 ) {
+LAST:
+  if ( dpm_base_ordtype == 1 ) {
     if ( pos1 < pos2 ) return 1;
     else if ( pos1 > pos2 ) return -1;
     else return (*cmpdl)(n,d1,d2);
@@ -2718,24 +2877,33 @@ int compdmm(int n,DMM m1,DMM m2)
 {
   int t;
 
-  if ( dpm_ordtype == 1 ) {
-    if ( m1->pos < m2->pos ) return 1;
-    else if ( m1->pos > m2->pos ) return -1;
-    else return (*cmpdl)(n,m1->dl,m2->dl);
-  } else if ( dpm_ordtype == 0 ) {
+  switch ( dpm_ordtype ) {
+  case 0: /* TOP ord->pos */
     t = (*cmpdl)(n,m1->dl,m2->dl);
     if ( t ) return t;
     else if ( m1->pos < m2->pos ) return 1;
     else if ( m1->pos > m2->pos ) return -1;
     else return 0;
-  } else if ( dpm_ordtype == 2 ) {
+  case 1: /* POT : pos->ord */
+    if ( m1->pos < m2->pos ) return 1;
+    else if ( m1->pos > m2->pos ) return -1;
+    else return (*cmpdl)(n,m1->dl,m2->dl);
+  case 2:  /* wPOT: weight->pos->ord */
+    if ( m1->dl->td > m2->dl->td ) return 1;
+    else if ( m1->dl->td < m2->dl->td ) return 1;
+    else if ( m1->pos < m2->pos ) return 1;
+    else if ( m1->pos > m2->pos ) return -1;
+    else return (*cmpdl)(n,m1->dl,m2->dl);
+  case 3: /* Schreyer */
     return compdmm_schreyer(n,m1,m2);
+  default:
+    error("compdmm : invalid dpm_ordtype");
   }
 }
 
 void adddpm(VL vl,DPM p1,DPM p2,DPM *pr)
 {
-  int n;
+  int n,s;
   DMM m1,m2,mr=0,mr0;
   Obj t;
   DL d;
@@ -2745,8 +2913,9 @@ void adddpm(VL vl,DPM p1,DPM p2,DPM *pr)
   else if ( !p2 )
     *pr = p1;
   else {
-    for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
-      switch ( compdmm(n,m1,m2) ) {
+    for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
+      s = compdmm(n,m1,m2);
+      switch ( s ) {
         case 0:
           arf_add(vl,C(m1),C(m2),&t);
           if ( t ) {
@@ -2760,6 +2929,7 @@ void adddpm(VL vl,DPM p1,DPM p2,DPM *pr)
           NEXTDMM(mr0,mr); mr->pos = m2->pos; mr->dl = m2->dl; C(mr) = C(m2);
           m2 = NEXT(m2); break;
       }
+    }
     if ( !mr0 )
       if ( m1 )
         mr0 = m1;
@@ -2809,6 +2979,26 @@ void chsgndpm(DPM p,DPM *pr)
   }
 }
 
+void mulcmp(Obj c,MP m)
+{
+  MP t;
+  Obj c1;
+
+  for ( t = m; t; t = NEXT(t) ) {
+    arf_mul(CO,c,C(t),&c1); C(t) = c1;
+  }
+}
+
+void mulcdmm(Obj c,DMM m)
+{
+  DMM t;
+  Obj c1;
+
+  for ( t = m; t; t = NEXT(t) ) {
+    arf_mul(CO,c,C(t),&c1); C(t) = c1;
+  }
+}
+
 void mulcdpm(VL vl,Obj c,DPM p,DPM *pr)
 {
   DMM m,mr=0,mr0;
@@ -2963,3 +3153,107 @@ int compdpm(VL vl,DPM p1,DPM p2)
   }
 }
 
+// p = ...+c*<<0,...0:pos>>+...
+DPM dpm_eliminate_term(DPM a,DPM p,Obj c,int pos)
+{
+  MP d0,d;
+  DMM m;
+  DP f;
+  DPM a1,p1,r;
+
+  if ( !a ) return 0;
+  d0 = 0;
+  for ( m = BDY(a); m; m = NEXT(m) )
+    if ( m->pos == pos ) {
+      NEXTMP(d0,d); d->dl = m->dl; arf_chsgn(m->c,&d->c);
+    }
+  if ( d0 ) {
+    NEXT(d) = 0; MKDP(NV(a),d0,f);
+    mulcdpm(CO,c,a,&a1);
+    mulobjdpm(CO,(Obj)f,p,&p1); 
+    adddpm(CO,a1,p1,&r);
+    return r;
+  } else
+    return a;
+}
+
+// <<...:i>> -> <<...:tab[i]>>
+DPM dpm_compress(DPM p,int *tab)
+{
+  DMM m,mr0,mr;
+  DPM t;
+
+  if ( !p ) return 0;
+  else {
+    for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
+      NEXTDMM(mr0,mr);
+      mr->dl = m->dl; mr->c = m->c; mr->pos = tab[m->pos];
+      if ( mr->pos <= 0 )
+        error("dpm_compress : invalid position");
+    }
+    NEXT(mr) = 0;
+    MKDPM(p->nv,mr0,t); t->sugar = p->sugar;
+    return t;
+  }
+}
+
+// input : s, s = syz(m) output simplified s, m
+void dpm_simplify_syz(LIST s,LIST m,LIST *s1,LIST *m1)
+{
+  int lm,ls,i,j,k,pos,nv;
+  DPM *am,*as;
+  DPM p;
+  DMM d;
+  Obj c;
+  int *tab,*dd;
+  NODE t,t1;
+
+  lm = length(BDY(m));
+  am = (DPM *)MALLOC((lm+1)*sizeof(DPM));
+  ls = length(BDY(s));
+  as = (DPM *)MALLOC(ls*sizeof(DPM));
+  for ( i = 1, t = BDY(m); i <= lm; i++, t = NEXT(t) ) am[i] = (DPM)BDY(t);
+  for ( i = 0, t = BDY(s); i < ls; i++, t = NEXT(t) ) as[i] = (DPM)BDY(t);
+
+  for ( i = 0; i < ls; i++ ) {
+    p = as[i];
+    if ( p == 0 ) continue;
+    nv = NV(p);
+    for ( d = BDY(p); d; d = NEXT(d) ) {
+      dd = d->dl->d;
+      for ( k = 0; k < nv; k++ ) if ( dd[k] ) break;
+      if ( k == nv ) break;
+    }
+    if ( d ) {
+      c = d->c; pos = d->pos;
+      for ( j = 0; j < ls; j++ )
+        if ( j != i ) {
+          as[j] = dpm_eliminate_term(as[j],p,c,pos);
+        }
+      // remove m[i]
+      am[pos] = 0;
+      as[i] = 0;
+    }
+  }
+  // compress s
+  // create index table from am[]
+  // (0 0 * 0 * ...) -> (0 0 1 0 2 ... ) which means 2->1, 4->2, ...
+  tab = (int *)MALLOC((lm+1)*sizeof(int));
+  for ( j = 0, i = 1; i <= lm; i++ ) {
+    if ( am[i] ) { j++; tab[i] = j; }
+    else tab[i] = 0;
+  }
+  t = 0;
+  for ( i = ls-1; i >= 0; i-- )
+    if ( as[i] ) {
+      p = dpm_compress(as[i],tab);
+      MKNODE(t1,(pointer)p,t); t = t1;
+    }
+  MKLIST(*s1,t);
+  t = 0;
+  for ( i = lm; i >= 1; i-- )
+    if ( am[i] ) {
+      MKNODE(t1,(pointer)am[i],t); t = t1;
+    }
+  MKLIST(*m1,t);
+}