[BACK]Return to nd.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2018 / engine

Diff for /OpenXM_contrib2/asir2018/engine/nd.c between version 1.21 and 1.38

version 1.21, 2019/09/19 06:29:48 version 1.38, 2020/10/26 02:41:05
Line 1 
Line 1 
 /* $OpenXM: OpenXM_contrib2/asir2018/engine/nd.c,v 1.20 2019/09/15 08:46:12 noro Exp $ */  /* $OpenXM: OpenXM_contrib2/asir2018/engine/nd.c,v 1.37 2020/10/06 06:31:19 noro Exp $ */
   
 #include "nd.h"  #include "nd.h"
   
Line 18  NM _nm_free_list;
Line 18  NM _nm_free_list;
 ND _nd_free_list;  ND _nd_free_list;
 ND_pairs _ndp_free_list;  ND_pairs _ndp_free_list;
 NODE nd_hcf;  NODE nd_hcf;
   int Nsyz,Nsamesig;
   
 Obj nd_top_weight;  Obj nd_top_weight;
   
Line 49  static NDV *nd_ps_trace;
Line 50  static NDV *nd_ps_trace;
 static NDV *nd_ps_sym;  static NDV *nd_ps_sym;
 static NDV *nd_ps_trace_sym;  static NDV *nd_ps_trace_sym;
 static RHist *nd_psh;  static RHist *nd_psh;
 static int nd_psn,nd_pslen;  static int nd_psn,nd_pslen,nd_nbase;
 static RHist *nd_red;  static RHist *nd_red;
 static int *nd_work_vector;  static int *nd_work_vector;
 static int **nd_matrix;  static int **nd_matrix;
Line 65  static int *nd_poly_weight,*nd_module_weight;
Line 66  static int *nd_poly_weight,*nd_module_weight;
 static NODE nd_tracelist;  static NODE nd_tracelist;
 static NODE nd_alltracelist;  static NODE nd_alltracelist;
 static int nd_gentrace,nd_gensyz,nd_nora,nd_newelim,nd_intersect,nd_lf;  static int nd_gentrace,nd_gensyz,nd_nora,nd_newelim,nd_intersect,nd_lf;
   static int nd_f4_td,nd_sba_f4step,nd_sba_pot,nd_sba_largelcm;
 static int *nd_gbblock;  static int *nd_gbblock;
 static NODE nd_nzlist,nd_check_splist;  static NODE nd_nzlist,nd_check_splist;
 static int nd_splist;  static int nd_splist;
 static int *nd_sugarweight;  static int *nd_sugarweight;
 static int nd_f4red,nd_rank0,nd_last_nonzero;  static int nd_f4red,nd_rank0,nd_last_nonzero;
   static DL *nd_sba_hm;
   static NODE *nd_sba_pos;
   
 NumberField get_numberfield();  NumberField get_numberfield();
 UINT *nd_det_compute_bound(NDV **dm,int n,int j);  UINT *nd_det_compute_bound(NDV **dm,int n,int j);
Line 744  int ndl_module_glex_compare(UINT *d1,UINT *d2)
Line 748  int ndl_module_glex_compare(UINT *d1,UINT *d2)
   
     default:      default:
       error("ndl_module_glex_compare : invalid module_ordtype");        error("ndl_module_glex_compare : invalid module_ordtype");
         return 0;
   }    }
 }  }
   
Line 754  int ndl_module_compare(UINT *d1,UINT *d2)
Line 759  int ndl_module_compare(UINT *d1,UINT *d2)
   
   switch ( nd_module_ordtype ) {    switch ( nd_module_ordtype ) {
     case 0:      case 0:
       if ( (c = ndl_lex_compare(d1,d2)) != 0 ) return c;        if ( (c = (*ndl_base_compare_function)(d1,d2)) != 0 ) return c;
       else if ( MPOS(d1) > MPOS(d2) ) return -1;        else if ( MPOS(d1) > MPOS(d2) ) return -1;
       else if ( MPOS(d1) < MPOS(d2) ) return 1;        else if ( MPOS(d1) < MPOS(d2) ) return 1;
       else return 0;        else return 0;
Line 776  int ndl_module_compare(UINT *d1,UINT *d2)
Line 781  int ndl_module_compare(UINT *d1,UINT *d2)
   
     default:      default:
       error("ndl_module_compare : invalid module_ordtype");        error("ndl_module_compare : invalid module_ordtype");
         return 0;
   }    }
 }  }
   
 extern DMMstack dmm_stack;  extern DMMstack dmm_stack;
 void _addtodl(int n,DL d1,DL d2);  void _addtodl(int n,DL d1,DL d2);
   void _adddl(int n,DL d1,DL d2,DL d3);
 int _eqdl(int n,DL d1,DL d2);  int _eqdl(int n,DL d1,DL d2);
   
 int ndl_module_schreyer_compare(UINT *m1,UINT *m2)  int ndl_module_schreyer_compare(UINT *m1,UINT *m2)
Line 840  LAST:
Line 847  LAST:
       break;        break;
     default:      default:
       error("ndl_schreyer_compare : invalid base ordtype");        error("ndl_schreyer_compare : invalid base ordtype");
         return 0;
   }    }
 }  }
   
Line 1197  INLINE int ndl_find_reducer(UINT *dg)
Line 1205  INLINE int ndl_find_reducer(UINT *dg)
     return -1;      return -1;
 }  }
   
   // ret=0,...,nd_psn-1 => reducer found
   // ret=nd_psn => reducer not found
   // ret=-1 => singular top reducible
   
   int comp_sig(SIG s1,SIG s2);
   void _ndltodl(UINT *ndl,DL dl);
   
   void print_sig(SIG s)
   {
     int i;
   
     fprintf(asir_out,"<<");
     for ( i = 0; i < nd_nvar; i++ ) {
       fprintf(asir_out,"%d",s->dl->d[i]);
       if ( i != nd_nvar-1 ) fprintf(asir_out,",");
     }
     fprintf(asir_out,">>*e%d",s->pos);
   }
   
   // assuming increasing order wrt signature
   
   INLINE int ndl_find_reducer_s(UINT *dg,SIG sig)
   {
     RHist r;
     int i,singular,ret,d,k;
     static int wpd,nvar;
     static SIG quo;
     static UINT *tmp;
   
     if ( !quo || nvar != nd_nvar ) NEWSIG(quo);
     if ( wpd != nd_wpd ) {
       wpd = nd_wpd;
       tmp = (UINT *)MALLOC(wpd*sizeof(UINT));
     }
     d = ndl_hash_value(dg);
   #if 1
     for ( r = nd_red[d], k = 0; r; r = NEXT(r), k++ ) {
       if ( ndl_equal(dg,DL(r)) ) {
         return r->index;
       }
     }
   #endif
     singular = 0;
     for ( i = 0; i < nd_psn; i++ ) {
       r = nd_psh[i];
       if ( ndl_reducible(dg,DL(r)) ) {
         ndl_sub(dg,DL(r),tmp);
         _ndltodl(tmp,DL(quo));
         _addtodl(nd_nvar,DL(nd_psh[i]->sig),DL(quo));
         quo->pos = nd_psh[i]->sig->pos;
         ret = comp_sig(sig,quo);
         if ( ret > 0 ) { singular = 0; break; }
         if ( ret == 0 ) { /* fprintf(asir_out,"s"); fflush(asir_out); */ singular = 1; }
       }
     }
     if ( singular ) return -1;
     else if ( i < nd_psn )
       nd_append_red(dg,i);
     return i;
   }
   
 ND nd_merge(ND p1,ND p2)  ND nd_merge(ND p1,ND p2)
 {  {
     int n,c;      int n,c;
Line 1539  int nd_nf(int mod,ND d,ND g,NDV *ps,int full,ND *rp)
Line 1608  int nd_nf(int mod,ND d,ND g,NDV *ps,int full,ND *rp)
     return 1;      return 1;
 }  }
   
   // ret=1 => success
   // ret=0 => overflow
   // ret=-1 => singular top reducible
   
   int nd_nf_s(int mod,ND d,ND g,NDV *ps,int full,ND *rp)
   {
       NM m,mrd,tail;
       NM mul;
       int n,sugar,psugar,sugar0,stat,index;
       int c,c1,c2,dummy;
       RHist h;
       NDV p,red;
       Q cg,cred,gcd,tq,qq;
       Z iq;
       DP dmul;
       NODE node;
       LIST hist;
       double hmag;
       P tp,tp1;
       Obj tr,tr1,div;
       union oNDC hg;
       P cont;
       SIG sig;
   
       if ( !g ) {
           *rp = d;
           return 1;
       }
       if ( !mod ) hmag = ((double)p_mag(HCP(g)))*nd_scale;
   
       sugar0 = sugar = SG(g);
       n = NV(g);
       mul = (NM)MALLOC(sizeof(struct oNM)+(nd_wpd-1)*sizeof(UINT));
       if ( d )
           for ( tail = BDY(d); NEXT(tail); tail = NEXT(tail) );
       sig = g->sig;
       for ( ; g; ) {
           index = ndl_find_reducer_s(HDL(g),sig);
           if ( index >= 0 && index < nd_psn ) {
               // reducer found
               h = nd_psh[index];
               ndl_sub(HDL(g),DL(h),DL(mul));
               if ( ndl_check_bound2(index,DL(mul)) ) {
                   nd_free(g); nd_free(d);
                   return 0;
               }
               p = ps[index];
               /* d+g -> div*(d+g)+mul*p */
               g = nd_reduce2(mod,d,g,p,mul,0,&div);
               sugar = MAX(sugar,SG(p)+TD(DL(mul)));
               if ( !mod && g && ((double)(p_mag(HCP(g))) > hmag) ) {
                   hg = HCU(g);
                   nd_removecont2(d,g);
                   hmag = ((double)p_mag(HCP(g)))*nd_scale;
               }
           } else if ( index == -1 ) {
             // singular top reducible
             return -1;
           } else if ( !full ) {
               *rp = g;
               g->sig = sig;
               return 1;
           } else {
               m = BDY(g);
               if ( NEXT(m) ) {
                   BDY(g) = NEXT(m); NEXT(m) = 0; LEN(g)--;
               } else {
                   FREEND(g); g = 0;
               }
               if ( d ) {
                   NEXT(tail)=m; tail=m; LEN(d)++;
               } else {
                   MKND(n,m,1,d); tail = BDY(d);
               }
           }
       }
       if ( d ) {
         SG(d) = sugar;
         d->sig = sig;
       }
       *rp = d;
       return 1;
   }
   
 int nd_nf_pbucket(int mod,ND g,NDV *ps,int full,ND *rp)  int nd_nf_pbucket(int mod,ND g,NDV *ps,int full,ND *rp)
 {  {
     int hindex,index;      int hindex,index;
Line 1651  int nd_nf_pbucket(int mod,ND g,NDV *ps,int full,ND *rp
Line 1804  int nd_nf_pbucket(int mod,ND g,NDV *ps,int full,ND *rp
     }      }
 }  }
   
   int nd_nf_pbucket_s(int mod,ND g,NDV *ps,int full,ND *rp)
   {
     int hindex,index;
     NDV p;
     ND u,d,red;
     NODE l;
     NM mul,m,mrd,tail;
     int sugar,psugar,n,h_reducible;
     PGeoBucket bucket;
     int c,c1,c2;
     Z cg,cred,gcd,zzz;
     RHist h;
     double hmag,gmag;
     int count = 0;
     int hcount = 0;
     SIG sig;
   
     if ( !g ) {
       *rp = 0;
       return 1;
     }
     sugar = SG(g);
     n = NV(g);
     if ( !mod ) hmag = ((double)p_mag((P)HCZ(g)))*nd_scale;
     bucket = create_pbucket();
     add_pbucket(mod,bucket,g);
     d = 0;
     mul = (NM)MALLOC(sizeof(struct oNM)+(nd_wpd-1)*sizeof(UINT));
     sig = g->sig;
     while ( 1 ) {
       if ( mod > 0 || mod == -1 )
         hindex = head_pbucket(mod,bucket);
       else if ( mod == -2 )
         hindex = head_pbucket_lf(bucket);
       else
         hindex = head_pbucket_q(bucket);
       if ( hindex < 0 ) {
         if ( DP_Print > 3 ) printf("(%d %d)",count,hcount);
         if ( d ) {
           SG(d) = sugar;
           d->sig = sig;
         }
         *rp = d;
         return 1;
       }
       g = bucket->body[hindex];
       index = ndl_find_reducer_s(HDL(g),sig);
       if ( index >= 0 && index < nd_psn ) {
         count++;
         if ( !d ) hcount++;
         h = nd_psh[index];
         ndl_sub(HDL(g),DL(h),DL(mul));
         if ( ndl_check_bound2(index,DL(mul)) ) {
           nd_free(d);
           free_pbucket(bucket);
           *rp = 0;
           return 0;
         }
         p = ps[index];
         if ( mod == -1 )
           CM(mul) = _mulsf(_invsf(HCM(p)),_chsgnsf(HCM(g)));
         else if ( mod == -2 ) {
           Z inv,t;
           divlf(ONE,HCZ(p),&inv);
           chsgnlf(HCZ(g),&t);
           mullf(inv,t,&CZ(mul));
         } else if ( mod ) {
           c1 = invm(HCM(p),mod); c2 = mod-HCM(g);
           DMAR(c1,c2,0,mod,c); CM(mul) = c;
         } else {
           igcd_cofactor(HCZ(g),HCZ(p),&gcd,&cg,&cred);
           chsgnz(cg,&CZ(mul));
           nd_mul_c_q(d,(P)cred);
           mulq_pbucket(bucket,cred);
           g = bucket->body[hindex];
           gmag = (double)p_mag((P)HCZ(g));
         }
         red = ndv_mul_nm(mod,mul,p);
         bucket->body[hindex] = nd_remove_head(g);
         red = nd_remove_head(red);
         add_pbucket(mod,bucket,red);
         psugar = SG(p)+TD(DL(mul));
         sugar = MAX(sugar,psugar);
         if ( !mod && hmag && (gmag > hmag) ) {
            g = normalize_pbucket(mod,bucket);
            if ( !g ) {
              if ( d ) {
                SG(d) = sugar;
                d->sig = sig;
              }
              *rp = d;
              return 1;
            }
            nd_removecont2(d,g);
            hmag = ((double)p_mag((P)HCZ(g)))*nd_scale;
            add_pbucket(mod,bucket,g);
         }
       } else if ( index == -1 ) {
         // singular top reducible
         return -1;
       } else if ( !full ) {
         g = normalize_pbucket(mod,bucket);
         if ( g ) {
           SG(g) = sugar;
           g->sig = sig;
         }
         *rp = g;
         return 1;
       } else {
         m = BDY(g);
         if ( NEXT(m) ) {
           BDY(g) = NEXT(m); NEXT(m) = 0; LEN(g)--;
         } else {
           FREEND(g); g = 0;
         }
         bucket->body[hindex] = g;
         NEXT(m) = 0;
         if ( d ) {
           NEXT(tail)=m; tail=m; LEN(d)++;
         } else {
           MKND(n,m,1,d); tail = BDY(d);
         }
       }
     }
   }
   
 /* input : list of NDV, cand : list of NDV */  /* input : list of NDV, cand : list of NDV */
   
 int ndv_check_membership(int m,NODE input,int obpe,int oadv,EPOS oepos,NODE cand)  int ndv_check_membership(int m,NODE input,int obpe,int oadv,EPOS oepos,NODE cand)
Line 1663  int ndv_check_membership(int m,NODE input,int obpe,int
Line 1942  int ndv_check_membership(int m,NODE input,int obpe,int
     Z q;      Z q;
     LIST list;      LIST list;
   
     ndv_setup(m,0,cand,nd_gentrace?1:0,1);      ndv_setup(m,0,cand,nd_gentrace?1:0,1,0);
     n = length(cand);      n = length(cand);
   
   if ( nd_gentrace ) { nd_alltracelist = 0; nd_tracelist = 0; }    if ( nd_gentrace ) { nd_alltracelist = 0; nd_tracelist = 0; }
Line 2111  NODE nd_gb(int m,int ishomo,int checkonly,int gensyz,i
Line 2390  NODE nd_gb(int m,int ishomo,int checkonly,int gensyz,i
   int diag_count = 0;    int diag_count = 0;
   P cont;    P cont;
   LIST list;    LIST list;
   struct oEGT eg1,eg2,eg_update;
   
   init_eg(&eg_update);
   Nnd_add = 0;    Nnd_add = 0;
   g = 0; d = 0;    g = 0; d = 0;
   for ( i = 0; i < nd_psn; i++ ) {    for ( i = 0; i < nd_psn; i++ ) {
Line 2176  again:
Line 2457  again:
         }          }
       }        }
       nfv = ndtondv(m,nf); nd_free(nf);        nfv = ndtondv(m,nf); nd_free(nf);
       nh = ndv_newps(m,nfv,0,0);        nh = ndv_newps(m,nfv,0);
       if ( !m && (ishomo && ++diag_count == diag_period) ) {        if ( !m && (ishomo && ++diag_count == diag_period) ) {
         diag_count = 0;          diag_count = 0;
         stat = do_diagonalize(sugar,m);          stat = do_diagonalize(sugar,m);
Line 2186  again:
Line 2467  again:
           goto again;            goto again;
         }          }
       }        }
   get_eg(&eg1);
       d = update_pairs(d,g,nh,0);        d = update_pairs(d,g,nh,0);
   get_eg(&eg2); add_eg(&eg_update,&eg1,&eg2);
       g = update_base(g,nh);        g = update_base(g,nh);
       FREENDP(l);        FREENDP(l);
    } else {     } else {
Line 2199  again:
Line 2482  again:
      }       }
      if ( DP_Print ) { printf("."); fflush(stdout); }       if ( DP_Print ) { printf("."); fflush(stdout); }
        FREENDP(l);         FREENDP(l);
        }
     }
     conv_ilist(nd_demand,0,g,indp);
     if ( !checkonly && DP_Print ) { printf("nd_gb done. Number of nd_add=%d\n",Nnd_add); fflush(stdout); }
   
     if ( DP_Print )
       print_eg("update",&eg_update);
     return g;
   }
   
   ND_pairs update_pairs_s(ND_pairs d,int t,NODE *syz);
   ND_pairs nd_newpairs_s(int t ,NODE *syz);
   
   int nd_nf_pbucket_s(int mod,ND g,NDV *ps,int full,ND *nf);
   int nd_nf_s(int mod,ND d,ND g,NDV *ps,int full,ND *nf);
   
   void _copydl(int n,DL d1,DL d2);
   void _subfromdl(int n,DL d1,DL d2);
   extern int (*cmpdl)(int n,DL d1,DL d2);
   
   NODE insert_sig(NODE l,SIG s)
   {
     int pos;
     DL sig;
     struct oNODE root;
     NODE p,prev,r;
     SIG t;
   
     pos = s->pos; sig = DL(s);
     root.next = l; prev = &root;
     for ( p = l; p; p = p->next ) {
       t = (SIG)p->body;
       if ( t->pos == pos ) {
         if ( _dl_redble(DL(t),sig,nd_nvar) )
           return root.next;
         else if ( _dl_redble(sig,DL(t),nd_nvar) )
           // remove p
           prev->next = p->next;
       } else
         prev = p;
     }
     NEWNODE(r); r->body = (pointer)s; r->next = 0;
     for ( p = &root; p->next; p = p->next );
     p->next = r;
   //  r->next = root.next;
   //  return r;
     return root.next;
   }
   
   ND_pairs remove_spair_s(ND_pairs d,SIG sig)
   {
     struct oND_pairs root;
     ND_pairs prev,p;
     SIG spsig;
   
     root.next = d;
     prev = &root; p = d;
     while ( p ) {
       spsig = p->sig;
       if ( sig->pos == spsig->pos && _dl_redble(DL(sig),DL(spsig),nd_nvar) ) {
         // remove p
         prev->next = p->next;
         Nsyz++;
       } else
         prev = p;
       p = p->next;
     }
     return (ND_pairs)root.next;
   }
   
   int _dl_redble_ext(DL,DL,DL,int);
   
   int small_lcm(ND_pairs l)
   {
     SIG sig;
     int i;
     NODE t;
     static DL lcm,mul,quo;
     static int nvar;
   
     if ( nd_sba_largelcm ) return 0;
     if ( nvar < nd_nvar ) {
       nvar = nd_nvar; NEWDL(lcm,nvar); NEWDL(quo,nvar); NEWDL(mul,nvar);
     }
     sig = l->sig;
     _ndltodl(l->lcm,lcm);
   #if 0
     for ( i = 0; i < nd_psn; i++ ) {
       if ( sig->pos == nd_psh[i]->sig->pos &&
         _dl_redble_ext(DL(nd_psh[i]->sig),DL(sig),quo,nd_nvar) ) {
         _ndltodl(DL(nd_psh[i]),mul);
         _addtodl(nd_nvar,quo,mul);
         if ( (*cmpdl)(nd_nvar,lcm,mul) > 0 )
           break;
       }
     }
     if ( i < nd_psn ) return 1;
     else return 0;
   #else
     for ( t = nd_sba_pos[sig->pos]; t; t = t->next ) {
       i = (long)BDY(t);
       if ( _dl_redble_ext(DL(nd_psh[i]->sig),DL(sig),quo,nd_nvar) ) {
         _ndltodl(DL(nd_psh[i]),mul);
         _addtodl(nd_nvar,quo,mul);
         if ( (*cmpdl)(nd_nvar,lcm,mul) > 0 )
           break;
       }
     }
     if ( t ) return 1;
     else return 0;
   #endif
   }
   
   ND_pairs remove_large_lcm(ND_pairs d)
   {
     struct oND_pairs root;
     ND_pairs prev,p;
   
     root.next = d;
     prev = &root; p = d;
     while ( p ) {
       if ( small_lcm(p) ) {
         // remove p
         prev->next = p->next;
       } else
         prev = p;
       p = p->next;
     }
     return (ND_pairs)root.next;
   }
   
   struct oEGT eg_create,eg_newpairs,eg_merge;
   
   NODE conv_ilist_s(int demand,int trace,int **indp);
   
   NODE nd_sba_buch(int m,int ishomo,int **indp)
   {
     int i,j,nh,sugar,stat;
     NODE r,t,g;
     ND_pairs d;
     ND_pairs l;
     ND h,nf,s,head,nf1;
     NDV nfv;
     Z q;
     union oNDC dn,hc;
     P cont;
     LIST list;
     SIG sig;
     NODE *syzlist;
     int Nredundant;
     DL lcm,quo,mul;
     struct oEGT eg1,eg2,eg_update,eg_remove,eg_large,eg_nf,eg_nfzero;
   
   init_eg(&eg_remove);
     syzlist = (NODE *)MALLOC(nd_psn*sizeof(NODE));
     Nsyz = 0;
     Nnd_add = 0;
     Nredundant = 0;
     d = 0;
     for ( i = 0; i < nd_psn; i++ ) {
       d = update_pairs_s(d,i,syzlist);
     }
     for ( i = 0; i < nd_psn; i++ )
       for ( j = i+1; j < nd_psn; j++ ) {
         NEWSIG(sig); sig->pos = j;
         _copydl(nd_nvar,nd_sba_hm[i],sig->dl);
         syzlist[sig->pos] = insert_sig(syzlist[sig->pos],sig);
       }
     sugar = 0;
     NEWDL(lcm,nd_nvar); NEWDL(quo,nd_nvar); NEWDL(mul,nd_nvar);
   init_eg(&eg_create);
   init_eg(&eg_merge);
   init_eg(&eg_large);
   init_eg(&eg_nf);
   init_eg(&eg_nfzero);
     while ( d ) {
   again:
       if ( DP_Print ) {
         int len;
         ND_pairs td;
         for ( td = d, len=0; td; td = td->next, len++)
           ;
          if ( !(len%100) ) fprintf(asir_out,"(%d)",len);
         }
       l = d; d = d->next;
       if ( small_lcm(l) ) {
         if ( DP_Print ) fprintf(asir_out,"M");
         Nredundant++;
         continue;
       }
       if ( SG(l) != sugar ) {
         sugar = SG(l);
         if ( DP_Print ) fprintf(asir_out,"%d",sugar);
       }
       sig = l->sig;
       stat = nd_sp(m,0,l,&h);
       if ( !stat ) {
         NEXT(l) = d; d = l;
         d = nd_reconstruct(0,d);
         goto again;
       }
   get_eg(&eg1);
   #if USE_GEOBUCKET
       stat = m?nd_nf_pbucket_s(m,h,nd_ps,!Top,&nf):nd_nf_s(m,0,h,nd_ps,!Top,&nf);
   #else
       stat = nd_nf_s(m,0,h,nd_ps,!Top,&nf);
   #endif
   get_eg(&eg2);
       if ( !stat ) {
         NEXT(l) = d; d = l;
         d = nd_reconstruct(0,d);
         goto again;
       } else if ( stat == -1 ) {
         if ( DP_Print ) { printf("S"); fflush(stdout); }
         FREENDP(l);
       } else if ( nf ) {
         if ( DP_Print ) { printf("+"); fflush(stdout); }
         add_eg(&eg_nf,&eg1,&eg2);
         hc = HCU(nf);
         nd_removecont(m,nf);
         nfv = ndtondv(m,nf); nd_free(nf);
         nh = ndv_newps(m,nfv,0);
   
         d = update_pairs_s(d,nh,syzlist);
         nd_sba_pos[sig->pos] = append_one(nd_sba_pos[sig->pos],nh);
         FREENDP(l);
      } else {
         add_eg(&eg_nfzero,&eg1,&eg2);
        // syzygy
   get_eg(&eg1);
        d = remove_spair_s(d,sig);
   get_eg(&eg2); add_eg(&eg_remove,&eg1,&eg2);
        syzlist[sig->pos] = insert_sig(syzlist[sig->pos],sig);
        if ( DP_Print ) { printf("."); fflush(stdout); }
        FREENDP(l);
    }     }
  }   }
  conv_ilist(nd_demand,0,g,indp);   g = conv_ilist_s(nd_demand,0,indp);
     if ( !checkonly && DP_Print ) { printf("nd_gb done. Number of nd_add=%d\n",Nnd_add); fflush(stdout); }   if ( DP_Print ) {
     return g;     printf("\nnd_sba done. nd_add=%d,Nsyz=%d,Nsamesig=%d,Nredundant=%d\n",Nnd_add,Nsyz,Nsamesig,Nredundant);
      fflush(stdout);
      print_eg("create",&eg_create);
      print_eg("merge",&eg_merge);
      print_eg("remove",&eg_remove);
      print_eg("nf",&eg_nf);
      print_eg("nfzero",&eg_nfzero);
      printf("\n");
    }
    return g;
 }  }
   
 /* splist = [[i1,i2],...] */  /* splist = [[i1,i2],...] */
Line 2459  again:
Line 2986  again:
             nd_tracelist = t;              nd_tracelist = t;
           }            }
         }          }
         nh = ndv_newps(0,nfv,nfqv,0);          nh = ndv_newps(0,nfv,nfqv);
         if ( ishomo && ++diag_count == diag_period ) {          if ( ishomo && ++diag_count == diag_period ) {
           diag_count = 0;            diag_count = 0;
           if ( DP_Print > 2 ) fprintf(asir_out,"|");            if ( DP_Print > 2 ) fprintf(asir_out,"|");
Line 2529  NODE ndv_reduceall(int m,NODE f)
Line 3056  NODE ndv_reduceall(int m,NODE f)
   
   if ( nd_nora ) return f;    if ( nd_nora ) return f;
   n = length(f);    n = length(f);
   ndv_setup(m,0,f,0,1);    ndv_setup(m,0,f,0,1,0);
   perm = (int *)MALLOC(n*sizeof(int));    perm = (int *)MALLOC(n*sizeof(int));
   if ( nd_gentrace ) {    if ( nd_gentrace ) {
     for ( t = nd_tracelist, i = 0; i < n; i++, t = NEXT(t) )      for ( t = nd_tracelist, i = 0; i < n; i++, t = NEXT(t) )
Line 2623  ND_pairs update_pairs( ND_pairs d, NODE /* of index */
Line 3150  ND_pairs update_pairs( ND_pairs d, NODE /* of index */
   }    }
 }  }
   
   ND_pairs merge_pairs_s(ND_pairs d,ND_pairs d1);
   
   ND_pairs update_pairs_s( ND_pairs d, int t,NODE *syz)
   {
     ND_pairs d1;
     struct oEGT eg1,eg2,eg3;
   
     if ( !t ) return d;
   get_eg(&eg1);
     d1 = nd_newpairs_s(t,syz);
   get_eg(&eg2); add_eg(&eg_create,&eg1,&eg2);
     d = merge_pairs_s(d,d1);
   get_eg(&eg3); add_eg(&eg_merge,&eg2,&eg3);
     return d;
   }
   
 ND_pairs nd_newpairs( NODE g, int t )  ND_pairs nd_newpairs( NODE g, int t )
 {  {
   NODE h;    NODE h;
Line 2661  ND_pairs nd_newpairs( NODE g, int t )
Line 3203  ND_pairs nd_newpairs( NODE g, int t )
   return r0;    return r0;
 }  }
   
   int comp_sig(SIG s1,SIG s2)
   {
     if ( nd_sba_pot ) {
       if ( s1->pos > s2->pos ) return 1;
       else if ( s1->pos < s2->pos ) return -1;
       else return (*cmpdl)(nd_nvar,s1->dl,s2->dl);
     } else {
       static DL m1,m2;
       static int nvar;
       int ret;
   
       if ( nvar != nd_nvar ) {
         nvar = nd_nvar; NEWDL(m1,nvar); NEWDL(m2,nvar);
       }
       _adddl(nd_nvar,s1->dl,nd_sba_hm[s1->pos],m1);
       _adddl(nd_nvar,s2->dl,nd_sba_hm[s2->pos],m2);
       ret = (*cmpdl)(nd_nvar,m1,m2);
       if ( ret != 0 ) return ret;
       else if ( s1->pos > s2->pos ) return 1;
       else if ( s1->pos < s2->pos ) return -1;
       else return 0;
     }
   }
   
   int _create_spair_s(int i1,int i2,ND_pairs sp,SIG sig1,SIG sig2)
   {
     int ret,s1,s2;
     RHist p1,p2;
     static int wpd;
     static UINT *lcm;
   
     sp->i1 = i1;
     sp->i2 = i2;
     p1 = nd_psh[i1];
     p2 = nd_psh[i2];
     ndl_lcm(DL(p1),DL(p2),sp->lcm);
     s1 = SG(p1)-TD(DL(p1));
     s2 = SG(p2)-TD(DL(p2));
     SG(sp) = MAX(s1,s2) + TD(sp->lcm);
   
     if ( wpd != nd_wpd ) {
       wpd = nd_wpd;
       lcm = (UINT *)MALLOC(wpd*sizeof(UINT));
     }
     // DL(sig1) <- sp->lcm
     // DL(sig1) -= DL(p1)
     // DL(sig1) += DL(p1->sig)
     ndl_sub(sp->lcm,DL(p1),lcm);
     _ndltodl(lcm,DL(sig1));
     _addtodl(nd_nvar,DL(p1->sig),DL(sig1));
     sig1->pos = p1->sig->pos;
   
     // DL(sig2) <- sp->lcm
     // DL(sig2) -= DL(p2)
     // DL(sig2) += DL(p2->sig)
     ndl_sub(sp->lcm,DL(p2),lcm);
     _ndltodl(lcm,DL(sig2));
     _addtodl(nd_nvar,DL(p2->sig),DL(sig2));
     sig2->pos = p2->sig->pos;
   
     ret = comp_sig(sig1,sig2);
     if ( ret == 0 ) return 0;
     else if ( ret > 0 ) sp->sig = sig1;
     else sp->sig = sig2;
     return 1;
   }
   
   SIG dup_sig(SIG sig)
   {
     SIG r;
   
     if ( !sig ) return 0;
     else {
       NEWSIG(r);
       _copydl(nd_nvar,DL(sig),DL(r));
       r->pos = sig->pos;
       return r;
     }
   }
   
   void dup_ND_pairs(ND_pairs to,ND_pairs from)
   {
     to->i1 = from->i1;
     to->i2 = from->i2;
     to->sugar = from->sugar;
     to->sugar2 = from->sugar2;
     ndl_copy(from->lcm,to->lcm);
     to->sig = dup_sig(from->sig);
   }
   
   ND_pairs merge_pairs_s(ND_pairs p1,ND_pairs p2)
   {
     struct oND_pairs root;
     ND_pairs q1,q2,r0,r;
     int ret;
   
     r = &root;
     for ( q1 = p1, q2 = p2; q1 != 0 && q2 != 0; ) {
       ret = comp_sig(q1->sig,q2->sig);
       if ( ret < 0 ) {
         r->next = q1; r = q1; q1 = q1->next;
       } else if ( ret > 0 ) {
         r->next = q2; r = q2; q2 = q2->next;
       } else {
         ret = DL_COMPARE(q1->lcm,q2->lcm);
         Nsamesig++;
         if ( ret < 0 ) {
           r->next = q1; r = q1; q1 = q1->next;
           q2 = q2->next;
         } else {
           r->next = q2; r = q2; q2 = q2->next;
           q1 = q1->next;
         }
       }
     }
     if ( q1 ) {
       r->next = q1;
     } else {
       r->next = q2;
     }
     return root.next;
   }
   
   ND_pairs insert_pair_s(ND_pairs l,ND_pairs s)
   {
     ND_pairs p,prev;
     int ret;
   
     for ( p = l, prev = 0; p != 0; prev = p, p = p->next ) {
       if ( (ret = comp_sig(s->sig,p->sig)) <= 0 )
         break;
     }
     if ( ret == 0 ) {
       ret = DL_COMPARE(s->lcm,p->lcm);
       if ( ret < 0 ) {
         // replace p with s
         s->next = p->next;
         if ( prev == 0 ) {
           return s;
         } else {
           prev->next = s;
           return l;
         }
       } else
         return l;
     } else {
       // insert s between prev and p
       s->next = p;
       if ( prev == 0 ) {
         return s;
       } else {
         prev->next = s;
         return l;
       }
     }
   }
   
   INLINE int __dl_redble(DL d1,DL d2,int nvar)
   {
     int i;
   
     if ( d1->td > d2->td )
       return 0;
     for ( i = nvar-1; i >= 0; i-- )
       if ( d1->d[i] > d2->d[i] )
         break;
     if ( i >= 0 )
       return 0;
     else
       return 1;
   }
   
   ND_pairs nd_newpairs_s(int t, NODE *syz)
   {
     NODE h,s;
     UINT *dl;
     int ts,ret,i;
     ND_pairs r,r0,_sp,sp;
     SIG spsig,tsig;
     static int nvar;
     static SIG _sig1,_sig2;
     struct oEGT eg1,eg2,eg3,eg4;
   
     NEWND_pairs(_sp);
     if ( !_sig1 || nvar != nd_nvar ) {
       nvar = nd_nvar; NEWSIG(_sig1); NEWSIG(_sig2);
     }
     r0 = 0;
     for ( i = 0; i < t; i++ ) {
       ret = _create_spair_s(i,t,_sp,_sig1,_sig2);
       if ( ret ) {
         spsig = _sp->sig;
         for ( s = syz[spsig->pos]; s; s = s->next ) {
           tsig = (SIG)s->body;
           if ( _dl_redble(DL(tsig),DL(spsig),nd_nvar) )
             break;
         }
         if ( s == 0 ) {
           NEWND_pairs(sp);
           dup_ND_pairs(sp,_sp);
           r0 = insert_pair_s(r0,sp);
         } else
           Nsyz++;
       }
     }
     return r0;
   }
   
 /* ipair = [i1,i2],[i1,i2],... */  /* ipair = [i1,i2],[i1,i2],... */
 ND_pairs nd_ipairtospair(NODE ipair)  ND_pairs nd_ipairtospair(NODE ipair)
 {  {
Line 2938  ND_pairs nd_minsugarp( ND_pairs d, ND_pairs *prest )
Line 3688  ND_pairs nd_minsugarp( ND_pairs d, ND_pairs *prest )
     return dm0;      return dm0;
 }  }
   
   ND_pairs nd_minsugarp_s( ND_pairs d, ND_pairs *prest )
   {
     int msugar;
     ND_pairs t,last;
   
     for ( msugar = SG(d), t = d; t; t = NEXT(t) )
       if ( SG(t) == msugar ) last = t;
     *prest = last->next;
     last->next = 0;
     return d;
   }
   
 int nd_tdeg(NDV c)  int nd_tdeg(NDV c)
 {  {
   int wmax = 0;    int wmax = 0;
Line 2950  int nd_tdeg(NDV c)
Line 3712  int nd_tdeg(NDV c)
   return wmax;    return wmax;
 }  }
   
 int ndv_newps(int m,NDV a,NDV aq,int f4)  int ndv_newps(int m,NDV a,NDV aq)
 {  {
     int len;      int len;
     RHist r;      RHist r;
Line 2986  int ndv_newps(int m,NDV a,NDV aq,int f4)
Line 3748  int ndv_newps(int m,NDV a,NDV aq,int f4)
         SG(r) = nd_tdeg(aq);          SG(r) = nd_tdeg(aq);
 #endif  #endif
         ndl_copy(HDL(aq),DL(r));          ndl_copy(HDL(aq),DL(r));
           r->sig = dup_sig(aq->sig);
     } else {      } else {
         if ( !m ) register_hcf(a);          if ( !m ) register_hcf(a);
         nd_bound[nd_psn] = ndv_compute_bound(a);          nd_bound[nd_psn] = ndv_compute_bound(a);
Line 2995  int ndv_newps(int m,NDV a,NDV aq,int f4)
Line 3758  int ndv_newps(int m,NDV a,NDV aq,int f4)
         SG(r) = nd_tdeg(a);          SG(r) = nd_tdeg(a);
 #endif  #endif
         ndl_copy(HDL(a),DL(r));          ndl_copy(HDL(a),DL(r));
           r->sig = dup_sig(a->sig);
     }      }
     if ( nd_demand ) {      if ( nd_demand ) {
         if ( aq ) {          if ( aq ) {
Line 3020  int ndv_newps(int m,NDV a,NDV aq,int f4)
Line 3784  int ndv_newps(int m,NDV a,NDV aq,int f4)
 /* nd_tracelist = [[0,index,div],...,[nd_psn-1,index,div]] */  /* nd_tracelist = [[0,index,div],...,[nd_psn-1,index,div]] */
 /* return 1 if success, 0 if failure (HC(a mod p)) */  /* return 1 if success, 0 if failure (HC(a mod p)) */
   
 int ndv_setup(int mod,int trace,NODE f,int dont_sort,int dont_removecont)  int ndv_setup(int mod,int trace,NODE f,int dont_sort,int dont_removecont,int sba)
 {  {
   int i,j,td,len,max;    int i,j,td,len,max;
   NODE s,s0,f0,tn;    NODE s,s0,f0,tn;
Line 3044  int ndv_setup(int mod,int trace,NODE f,int dont_sort,i
Line 3808  int ndv_setup(int mod,int trace,NODE f,int dont_sort,i
     if ( BDY(s) ) { w[i].p = BDY(s); w[i].i = j; i++; }      if ( BDY(s) ) { w[i].p = BDY(s); w[i].i = j; i++; }
   if ( !dont_sort ) {    if ( !dont_sort ) {
     /* XXX heuristic */      /* XXX heuristic */
     if ( !nd_ord->id && (nd_ord->ord.simple<2) )      if ( !sba && !nd_ord->id && (nd_ord->ord.simple<2) )
       qsort(w,nd_psn,sizeof(struct oNDVI),        qsort(w,nd_psn,sizeof(struct oNDVI),
         (int (*)(const void *,const void *))ndvi_compare_rev);          (int (*)(const void *,const void *))ndvi_compare_rev);
     else      else
Line 3120  int ndv_setup(int mod,int trace,NODE f,int dont_sort,i
Line 3884  int ndv_setup(int mod,int trace,NODE f,int dont_sort,i
       }        }
     }      }
   }    }
     if ( sba ) {
       nd_sba_hm = (DL *)MALLOC(nd_psn*sizeof(DL));
      // setup signatures
       for ( i = 0; i < nd_psn; i++ ) {
         SIG sig;
   
         NEWSIG(sig); sig->pos = i;
         nd_ps[i]->sig = sig;
         if ( nd_demand ) nd_ps_sym[i]->sig = sig;
           nd_psh[i]->sig = sig;
         if ( trace ) {
           nd_ps_trace[i]->sig = sig;
           if ( nd_demand ) nd_ps_trace_sym[i]->sig = sig;
         }
         NEWDL(nd_sba_hm[i],nd_nvar);
         _ndltodl(DL(nd_psh[i]),nd_sba_hm[i]);
       }
       nd_sba_pos = (NODE *)MALLOC(nd_psn*sizeof(NODE));
       for ( i = 0; i < nd_psn; i++ ) {
         j = nd_psh[i]->sig->pos;
         nd_sba_pos[j] = append_one(nd_sba_pos[j],i);
       }
     }
   if ( nd_gentrace && nd_tracelist ) NEXT(tn) = 0;    if ( nd_gentrace && nd_tracelist ) NEXT(tn) = 0;
   return 1;    return 1;
 }  }
Line 3356  void nd_gr(LIST f,LIST v,int m,int homo,int retdp,int 
Line 4143  void nd_gr(LIST f,LIST v,int m,int homo,int retdp,int 
             ndv_homogenize((NDV)BDY(t),obpe,oadv,oepos,ompos);              ndv_homogenize((NDV)BDY(t),obpe,oadv,oepos,ompos);
     }      }
   
     ndv_setup(m,0,fd0,(nd_gbblock||nd_splist||nd_check_splist)?1:0,0);      ndv_setup(m,0,fd0,(nd_gbblock||nd_splist||nd_check_splist)?1:0,0,0);
     if ( nd_gentrace ) {      if ( nd_gentrace ) {
         MKLIST(l1,nd_tracelist); MKNODE(nd_alltracelist,l1,0);          MKLIST(l1,nd_tracelist); MKNODE(nd_alltracelist,l1,0);
     }      }
Line 3459  FINAL:
Line 4246  FINAL:
 #endif  #endif
 }  }
   
   NODE nd_sba_f4(int m,int **indp);
   
   void nd_sba(LIST f,LIST v,int m,int homo,int retdp,int f4,struct order_spec *ord,LIST *rp)
   {
     VL tv,fv,vv,vc,av;
     NODE fd,fd0,r,r0,t,x,s,xx;
     int e,max,nvar,i;
     NDV b;
     int ishomo,nalg,wmax,len;
     NMV a;
     P p,zp;
     Q dmy;
     struct order_spec *ord1;
     int j;
     int *perm;
     EPOS oepos;
     int obpe,oadv,ompos,cbpe;
     struct oEGT eg0,eg1,egconv;
   
     nd_module = 0;
     nd_demand = 0;
     parse_nd_option(current_option);
     Nsamesig = 0;
     if ( DP_Multiple )
       nd_scale = ((double)DP_Multiple)/(double)(Denominator?Denominator:1);
     get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc);
     if ( m && nd_vc )
       error("nd_sba : computation over Fp(X) is unsupported. Use dp_gr_mod_main().");
     for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ );
     switch ( ord->id ) {
       case 1:
         if ( ord->nv != nvar )
           error("nd_sba : invalid order specification");
           break;
         default:
           break;
     }
     nd_nalg = 0;
     nd_init_ord(ord);
     // for SIG comparison
     initd(ord);
     for ( t = BDY(f), max = 1; t; t = NEXT(t) ) {
       for ( tv = vv; tv; tv = NEXT(tv) ) {
         e = getdeg(tv->v,(P)BDY(t));
         max = MAX(e,max);
       }
     }
     nd_setup_parameters(nvar,max);
     obpe = nd_bpe; oadv = nmv_adv; oepos = nd_epos; ompos = nd_mpos;
     ishomo = 1;
     for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
       if ( !m ) ptozp((P)BDY(t),1,&dmy,&zp);
       else zp = (P)BDY(t);
       b = (pointer)ptondv(CO,vv,zp);
       if ( ishomo )
         ishomo = ishomo && ndv_ishomo(b);
       if ( m ) ndv_mod(m,b);
       if ( b ) { NEXTNODE(fd0,fd); BDY(fd) = (pointer)b; }
     }
     if ( fd0 ) NEXT(fd) = 0;
   
     if ( !ishomo && homo ) {
       for ( t = fd0, wmax = max; t; t = NEXT(t) ) {
         b = (NDV)BDY(t); len = LEN(b);
         for ( a = BDY(b), i = 0; i < len; i++, NMV_ADV(a) )
           wmax = MAX(TD(DL(a)),wmax);
         }
         homogenize_order(ord,nvar,&ord1);
         nd_init_ord(ord1);
         nd_setup_parameters(nvar+1,nd_nzlist?0:wmax);
         for ( t = fd0; t; t = NEXT(t) )
           ndv_homogenize((NDV)BDY(t),obpe,oadv,oepos,ompos);
     }
   
     ndv_setup(m,0,fd0,0,0,1);
     x = f4 ? nd_sba_f4(m,&perm) : nd_sba_buch(m,ishomo || homo,&perm);
     if ( !x ) {
       *rp = 0; return;
     }
     if ( !ishomo && homo ) {
          /* dehomogenization */
       for ( t = x; t; t = NEXT(t) ) ndv_dehomogenize((NDV)BDY(t),ord);
       nd_init_ord(ord);
       nd_setup_parameters(nvar,0);
     }
     nd_demand = 0;
     x = ndv_reducebase(x,perm);
     x = ndv_reduceall(m,x);
     nd_setup_parameters(nd_nvar,0);
     get_eg(&eg0);
     for ( r0 = 0, t = x; t; t = NEXT(t) ) {
       NEXTNODE(r0,r);
       if ( retdp ) BDY(r) = ndvtodp(m,BDY(t));
       else BDY(r) = ndvtop(m,CO,vv,BDY(t));
     }
     if ( r0 ) NEXT(r) = 0;
     MKLIST(*rp,r0);
     get_eg(&eg1); init_eg(&egconv); add_eg(&egconv,&eg0,&eg1);
     print_eg("conv",&egconv); fprintf(asir_out,"\n");
   }
   
 void nd_gr_postproc(LIST f,LIST v,int m,struct order_spec *ord,int do_check,LIST *rp)  void nd_gr_postproc(LIST f,LIST v,int m,struct order_spec *ord,int do_check,LIST *rp)
 {  {
     VL tv,fv,vv,vc,av;      VL tv,fv,vv,vc,av;
Line 3518  void nd_gr_postproc(LIST f,LIST v,int m,struct order_s
Line 4406  void nd_gr_postproc(LIST f,LIST v,int m,struct order_s
         if ( b ) { NEXTNODE(fd0,fd); BDY(fd) = (pointer)b; }          if ( b ) { NEXTNODE(fd0,fd); BDY(fd) = (pointer)b; }
     }      }
     if ( fd0 ) NEXT(fd) = 0;      if ( fd0 ) NEXT(fd) = 0;
     ndv_setup(m,0,fd0,0,1);      ndv_setup(m,0,fd0,0,1,0);
     for ( x = 0, i = 0; i < nd_psn; i++ )      for ( x = 0, i = 0; i < nd_psn; i++ )
         x = update_base(x,i);          x = update_base(x,i);
     if ( do_check ) {      if ( do_check ) {
Line 3839  void nd_gr_trace(LIST f,LIST v,int trace,int homo,int 
Line 4727  void nd_gr_trace(LIST f,LIST v,int trace,int homo,int 
     tl1 = tl2 = tl3 = tl4 = 0;      tl1 = tl2 = tl3 = tl4 = 0;
         if ( Demand )          if ( Demand )
             nd_demand = 1;              nd_demand = 1;
         ret = ndv_setup(m,1,fd0,nd_gbblock?1:0,0);          ret = ndv_setup(m,1,fd0,nd_gbblock?1:0,0,0);
         if ( nd_gentrace ) {          if ( nd_gentrace ) {
             MKLIST(l1,nd_tracelist); MKNODE(nd_alltracelist,l1,0);              MKLIST(l1,nd_tracelist); MKNODE(nd_alltracelist,l1,0);
         }          }
Line 3972  DL ndltodl(int n,UINT *ndl)
Line 4860  DL ndltodl(int n,UINT *ndl)
     int i,j,l,s,ord_l;      int i,j,l,s,ord_l;
     struct order_pair *op;      struct order_pair *op;
   
     NEWDL(dl,n);      NEWDL_NOINIT(dl,n);
     dl->td = TD(ndl);      dl->td = TD(ndl);
     d = dl->d;      d = dl->d;
     if ( nd_blockmask ) {      if ( nd_blockmask ) {
Line 3988  DL ndltodl(int n,UINT *ndl)
Line 4876  DL ndltodl(int n,UINT *ndl)
     return dl;      return dl;
 }  }
   
   void _ndltodl(UINT *ndl,DL dl)
   {
       int *d;
       int i,j,l,s,ord_l,n;
       struct order_pair *op;
   
       n = nd_nvar;
       dl->td = TD(ndl);
       d = dl->d;
       if ( nd_blockmask ) {
           l = nd_blockmask->n;
           op = nd_blockmask->order_pair;
           for ( j = 0, s = 0; j < l; j++ ) {
               ord_l = op[j].length;
               for ( i = 0; i < ord_l; i++, s++ ) d[s] = GET_EXP(ndl,s);
           }
       } else {
           for ( i = 0; i < n; i++ ) d[i] = GET_EXP(ndl,i);
       }
   }
   
 void nmtodp(int mod,NM m,DP *r)  void nmtodp(int mod,NM m,DP *r)
 {  {
     DP dp;      DP dp;
Line 4045  void nd_print(ND p)
Line 4954  void nd_print(ND p)
     else {      else {
         for ( m = BDY(p); m; m = NEXT(m) ) {          for ( m = BDY(p); m; m = NEXT(m) ) {
             if ( CM(m) & 0x80000000 ) printf("+@_%d*",IFTOF(CM(m)));              if ( CM(m) & 0x80000000 ) printf("+@_%d*",IFTOF(CM(m)));
             else printf("+%d*",CM(m));              else printf("+%ld*",CM(m));
             ndl_print(DL(m));              ndl_print(DL(m));
         }          }
         printf("\n");          printf("\n");
Line 4640  ND_pairs nd_reconstruct(int trace,ND_pairs d)
Line 5549  ND_pairs nd_reconstruct(int trace,ND_pairs d)
         NEXTND_pairs(s0,s);          NEXTND_pairs(s0,s);
         s->i1 = t->i1;          s->i1 = t->i1;
         s->i2 = t->i2;          s->i2 = t->i2;
           s->sig = t->sig;
         SG(s) = SG(t);          SG(s) = SG(t);
         ndl_reconstruct(LCM(t),LCM(s),obpe,oepos);          ndl_reconstruct(LCM(t),LCM(s),obpe,oepos);
     }      }
Line 4658  ND_pairs nd_reconstruct(int trace,ND_pairs d)
Line 5568  ND_pairs nd_reconstruct(int trace,ND_pairs d)
             h = ndl_hash_value(DL(mr));              h = ndl_hash_value(DL(mr));
             NEXT(mr) = nd_red[h];              NEXT(mr) = nd_red[h];
             nd_red[h] = mr;              nd_red[h] = mr;
               mr->sig = r->sig;
         }          }
     for ( i = 0; i < REDTAB_LEN; i++ ) old_red[i] = 0;      for ( i = 0; i < REDTAB_LEN; i++ ) old_red[i] = 0;
     old_red = 0;      old_red = 0;
     for ( i = 0; i < nd_psn; i++ ) {      for ( i = 0; i < nd_psn; i++ ) {
         NEWRHist(r); SG(r) = SG(nd_psh[i]);          NEWRHist(r); SG(r) = SG(nd_psh[i]);
         ndl_reconstruct(DL(nd_psh[i]),DL(r),obpe,oepos);          ndl_reconstruct(DL(nd_psh[i]),DL(r),obpe,oepos);
           r->sig = nd_psh[i]->sig;
         nd_psh[i] = r;          nd_psh[i] = r;
     }      }
     if ( s0 ) NEXT(s) = 0;      if ( s0 ) NEXT(s) = 0;
Line 4675  ND_pairs nd_reconstruct(int trace,ND_pairs d)
Line 5587  ND_pairs nd_reconstruct(int trace,ND_pairs d)
     return s0;      return s0;
 }  }
   
   void nd_reconstruct_s(int trace,ND_pairs *d)
   {
       int i,obpe,oadv,h;
       static NM prev_nm_free_list;
       static ND_pairs prev_ndp_free_list;
       RHist mr0,mr;
       RHist r;
       RHist *old_red;
       ND_pairs s0,s,t;
       EPOS oepos;
   
       obpe = nd_bpe;
       oadv = nmv_adv;
       oepos = nd_epos;
       if ( obpe < 2 ) nd_bpe = 2;
       else if ( obpe < 3 ) nd_bpe = 3;
       else if ( obpe < 4 ) nd_bpe = 4;
       else if ( obpe < 5 ) nd_bpe = 5;
       else if ( obpe < 6 ) nd_bpe = 6;
       else if ( obpe < 8 ) nd_bpe = 8;
       else if ( obpe < 10 ) nd_bpe = 10;
       else if ( obpe < 16 ) nd_bpe = 16;
       else if ( obpe < 32 ) nd_bpe = 32;
       else error("nd_reconstruct_s : exponent too large");
   
       nd_setup_parameters(nd_nvar,0);
       prev_nm_free_list = _nm_free_list;
       prev_ndp_free_list = _ndp_free_list;
       _nm_free_list = 0;
       _ndp_free_list = 0;
       for ( i = nd_psn-1; i >= 0; i-- ) {
           ndv_realloc(nd_ps[i],obpe,oadv,oepos);
           ndv_realloc(nd_ps_sym[i],obpe,oadv,oepos);
       }
       if ( trace )
           for ( i = nd_psn-1; i >= 0; i-- ) {
               ndv_realloc(nd_ps_trace[i],obpe,oadv,oepos);
               ndv_realloc(nd_ps_trace_sym[i],obpe,oadv,oepos);
           }
   
       for ( i = 0; i < nd_nbase; i++ ) {
         s0 = 0;
         for ( t = d[i]; t; t = NEXT(t) ) {
             NEXTND_pairs(s0,s);
             s->i1 = t->i1;
             s->i2 = t->i2;
             s->sig = t->sig;
             SG(s) = SG(t);
             ndl_reconstruct(LCM(t),LCM(s),obpe,oepos);
         }
         d[i] = s0;
       }
   
       old_red = (RHist *)MALLOC(REDTAB_LEN*sizeof(RHist));
       for ( i = 0; i < REDTAB_LEN; i++ ) {
           old_red[i] = nd_red[i];
           nd_red[i] = 0;
       }
       for ( i = 0; i < REDTAB_LEN; i++ )
           for ( r = old_red[i]; r; r = NEXT(r) ) {
               NEWRHist(mr);
               mr->index = r->index;
               SG(mr) = SG(r);
               ndl_reconstruct(DL(r),DL(mr),obpe,oepos);
               h = ndl_hash_value(DL(mr));
               NEXT(mr) = nd_red[h];
               nd_red[h] = mr;
               mr->sig = r->sig;
           }
       for ( i = 0; i < REDTAB_LEN; i++ ) old_red[i] = 0;
       old_red = 0;
       for ( i = 0; i < nd_psn; i++ ) {
           NEWRHist(r); SG(r) = SG(nd_psh[i]);
           ndl_reconstruct(DL(nd_psh[i]),DL(r),obpe,oepos);
           r->sig = nd_psh[i]->sig;
           nd_psh[i] = r;
       }
       if ( s0 ) NEXT(s) = 0;
       prev_nm_free_list = 0;
       prev_ndp_free_list = 0;
   #if 0
       GC_gcollect();
   #endif
   }
   
 void ndl_reconstruct(UINT *d,UINT *r,int obpe,EPOS oepos)  void ndl_reconstruct(UINT *d,UINT *r,int obpe,EPOS oepos)
 {  {
     int n,i,ei,oepw,omask0,j,s,ord_l,l;      int n,i,ei,oepw,omask0,j,s,ord_l,l;
Line 4783  int nd_sp(int mod,int trace,ND_pairs p,ND *rp)
Line 5780  int nd_sp(int mod,int trace,ND_pairs p,ND *rp)
         MKLIST(hist,node); MKNODE(node,hist,nd_tracelist);          MKLIST(hist,node); MKNODE(node,hist,nd_tracelist);
         nd_tracelist = node;          nd_tracelist = node;
     }      }
       if ( *rp )
         (*rp)->sig = p->sig;
     FREENM(m1); FREENM(m2);      FREENM(m1); FREENM(m2);
     return 1;      return 1;
 }  }
Line 5161  NDV ndv_dup_realloc(NDV p,int obpe,int oadv,EPOS oepos
Line 6160  NDV ndv_dup_realloc(NDV p,int obpe,int oadv,EPOS oepos
     }      }
     MKNDV(NV(p),mr0,len,r);      MKNDV(NV(p),mr0,len,r);
     SG(r) = SG(p);      SG(r) = SG(p);
       r->sig = p->sig;
     return r;      return r;
 }  }
   
Line 5493  NDV ndtondv(int mod,ND p)
Line 6493  NDV ndtondv(int mod,ND p)
     }      }
     MKNDV(NV(p),m0,len,d);      MKNDV(NV(p),m0,len,d);
     SG(d) = SG(p);      SG(d) = SG(p);
       d->sig = p->sig;
     return d;      return d;
 }  }
   
Line 5538  void dpm_sort(DPM p,DPM *rp)
Line 6539  void dpm_sort(DPM p,DPM *rp)
   
 int dpm_comp(DPM *a,DPM *b)  int dpm_comp(DPM *a,DPM *b)
 {  {
   return compdpm(CO,*a,*b);    return -compdpm(CO,*a,*b);
 }  }
   
 NODE dpm_sort_list(NODE l)  NODE dpm_sort_list(NODE l)
Line 5616  ND ndvtond(int mod,NDV p)
Line 6617  ND ndvtond(int mod,NDV p)
     NEXT(m) = 0;      NEXT(m) = 0;
     MKND(NV(p),m0,len,d);      MKND(NV(p),m0,len,d);
     SG(d) = SG(p);      SG(d) = SG(p);
       d->sig = p->sig;
     return d;      return d;
 }  }
   
Line 5694  void ndv_print(NDV p)
Line 6696  void ndv_print(NDV p)
         len = LEN(p);          len = LEN(p);
         for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) {          for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) {
             if ( CM(m) & 0x80000000 ) printf("+@_%d*",IFTOF(CM(m)));              if ( CM(m) & 0x80000000 ) printf("+@_%d*",IFTOF(CM(m)));
             else printf("+%d*",CM(m));              else printf("+%ld*",CM(m));
             ndl_print(DL(m));              ndl_print(DL(m));
         }          }
         printf("\n");          printf("\n");
Line 6027  void nd_nf_p(Obj f,LIST g,LIST v,int m,struct order_sp
Line 7029  void nd_nf_p(Obj f,LIST g,LIST v,int m,struct order_sp
     ndf = (pointer)ndvtond(m,ndvf);      ndf = (pointer)ndvtond(m,ndvf);
   
     /* dont sort, dont removecont */      /* dont sort, dont removecont */
     ndv_setup(m,0,in0,1,1);      ndv_setup(m,0,in0,1,1,0);
     nd_scale=2;      nd_scale=2;
     stat = nd_nf(m,0,ndf,nd_ps,1,&nf);      stat = nd_nf(m,0,ndf,nd_ps,1,&nf);
     if ( !stat )      if ( !stat )
Line 6887  int nd_symbolic_preproc(PGeoBucket bucket,int trace,UI
Line 7889  int nd_symbolic_preproc(PGeoBucket bucket,int trace,UI
             if ( ndl_check_bound2(index,DL(mul)) )              if ( ndl_check_bound2(index,DL(mul)) )
                 return 0;                  return 0;
             sugar = TD(DL(mul))+SG(ps[index]);              sugar = TD(DL(mul))+SG(ps[index]);
             MKNM_ind_pair(pair,mul,index,sugar);              MKNM_ind_pair(pair,mul,index,sugar,0);
             red = ndv_mul_nm_symbolic(mul,ps[index]);              red = ndv_mul_nm_symbolic(mul,ps[index]);
             add_pbucket_symbolic(bucket,nd_remove_head(red));              add_pbucket_symbolic(bucket,nd_remove_head(red));
             NEXTNODE(rp0,rp); BDY(rp) = (pointer)pair;              NEXTNODE(rp0,rp); BDY(rp) = (pointer)pair;
Line 6992  NODE nd_f4(int m,int checkonly,int **indp)
Line 7994  NODE nd_f4(int m,int checkonly,int **indp)
         if ( nflist ) nd_last_nonzero = f4red;          if ( nflist ) nd_last_nonzero = f4red;
         for ( r = nflist; r; r = NEXT(r) ) {          for ( r = nflist; r; r = NEXT(r) ) {
             nf = (NDV)BDY(r);              nf = (NDV)BDY(r);
               if ( nd_f4_td ) SG(nf) = nd_tdeg(nf);
             ndv_removecont(m,nf);              ndv_removecont(m,nf);
             if ( !m && nd_nalg ) {              if ( !m && nd_nalg ) {
                 ND nf1;                  ND nf1;
Line 7001  NODE nd_f4(int m,int checkonly,int **indp)
Line 8004  NODE nd_f4(int m,int checkonly,int **indp)
                 nd_removecont(m,nf1);                  nd_removecont(m,nf1);
                 nf = ndtondv(m,nf1);                  nf = ndtondv(m,nf1);
             }              }
             nh = ndv_newps(m,nf,0,1);              nh = ndv_newps(m,nf,0);
             d = update_pairs(d,g,nh,0);              d = update_pairs(d,g,nh,0);
             g = update_base(g,nh);              g = update_base(g,nh);
         }          }
Line 7129  NODE nd_f4_trace(int m,int **indp)
Line 8132  NODE nd_f4_trace(int m,int **indp)
             nfv = ndv_dup(0,nfqv);              nfv = ndv_dup(0,nfqv);
             ndv_mod(m,nfv);              ndv_mod(m,nfv);
             ndv_removecont(m,nfv);              ndv_removecont(m,nfv);
             nh = ndv_newps(0,nfv,nfqv,1);              nh = ndv_newps(0,nfv,nfqv);
             d = update_pairs(d,g,nh,0);              d = update_pairs(d,g,nh,0);
             g = update_base(g,nh);              g = update_base(g,nh);
         }          }
Line 8540  void conv_ilist(int demand,int trace,NODE g,int **indp
Line 9543  void conv_ilist(int demand,int trace,NODE g,int **indp
   if ( indp ) *indp = ind;    if ( indp ) *indp = ind;
 }  }
   
   NODE conv_ilist_s(int demand,int trace,int **indp)
   {
     int n,i,j;
     int *ind;
     NODE g0,g;
   
     n = nd_psn;
     ind = (int *)MALLOC(n*sizeof(int));
     g0 = 0;
     for ( i = 0; i < n; i++ ) {
       ind[i] = i;
       NEXTNODE(g0,g);
       BDY(g) = (pointer)(demand?ndv_load(i):(trace?nd_ps_trace[i]:nd_ps[i]));
     }
     if ( g0 ) NEXT(g) = 0;
     if ( indp ) *indp = ind;
     return g0;
   }
   
 void parse_nd_option(NODE opt)  void parse_nd_option(NODE opt)
 {  {
     NODE t,p,u;    NODE t,p,u;
   int i,s,n;    int i,s,n;
     char *key;    char *key;
     Obj value;    Obj value;
   
     nd_gentrace = 0; nd_gensyz = 0; nd_nora = 0; nd_gbblock = 0;    nd_gentrace = 0; nd_gensyz = 0; nd_nora = 0; nd_gbblock = 0;
   nd_newelim = 0; nd_intersect = 0; nd_nzlist = 0;    nd_newelim = 0; nd_intersect = 0; nd_nzlist = 0;
   nd_splist = 0; nd_check_splist = 0;    nd_splist = 0; nd_check_splist = 0;
     nd_sugarweight = 0;    nd_sugarweight = 0; nd_f4red =0; nd_rank0 = 0;
     nd_f4red =0;    nd_f4_td = 0; nd_sba_f4step = 2; nd_sba_pot = 0; nd_sba_largelcm = 0;
     nd_rank0 = 0;    for ( t = opt; t; t = NEXT(t) ) {
     for ( t = opt; t; t = NEXT(t) ) {      p = BDY((LIST)BDY(t));
         p = BDY((LIST)BDY(t));      key = BDY((STRING)BDY(p));
         key = BDY((STRING)BDY(p));      value = (Obj)BDY(NEXT(p));
         value = (Obj)BDY(NEXT(p));      if ( !strcmp(key,"gentrace") )
         if ( !strcmp(key,"gentrace") )        nd_gentrace = value?1:0;
             nd_gentrace = value?1:0;      else if ( !strcmp(key,"gensyz") )
         else if ( !strcmp(key,"gensyz") )        nd_gensyz = value?1:0;
             nd_gensyz = value?1:0;      else if ( !strcmp(key,"nora") )
         else if ( !strcmp(key,"nora") )        nd_nora = value?1:0;
             nd_nora = value?1:0;      else if ( !strcmp(key,"gbblock") ) {
         else if ( !strcmp(key,"gbblock") ) {        if ( value && OID(value) == O_LIST ) {
             if ( value && OID(value) == O_LIST ) {  
         u = BDY((LIST)value);          u = BDY((LIST)value);
               nd_gbblock = MALLOC((2*length(u)+1)*sizeof(int));          nd_gbblock = MALLOC((2*length(u)+1)*sizeof(int));
         for ( i = 0; u; u = NEXT(u) ) {          for ( i = 0; u; u = NEXT(u) ) {
           p = BDY((LIST)BDY(u));            p = BDY((LIST)BDY(u));
           s = nd_gbblock[i++] = ZTOS((Q)BDY(p));            s = nd_gbblock[i++] = ZTOS((Q)BDY(p));
           nd_gbblock[i++] = s+ZTOS((Q)BDY(NEXT(p)))-1;            nd_gbblock[i++] = s+ZTOS((Q)BDY(NEXT(p)))-1;
         }          }
         nd_gbblock[i] = -1;          nd_gbblock[i] = -1;
             } else        } else
               nd_gbblock = 0;          nd_gbblock = 0;
     } else if ( !strcmp(key,"newelim") )      } else if ( !strcmp(key,"newelim") )
             nd_newelim = value?1:0;              nd_newelim = value?1:0;
     else if ( !strcmp(key,"intersect") )      else if ( !strcmp(key,"intersect") )
Line 8584  void parse_nd_option(NODE opt)
Line 9605  void parse_nd_option(NODE opt)
     else if ( !strcmp(key,"lf") )      else if ( !strcmp(key,"lf") )
             nd_lf = value?1:0;              nd_lf = value?1:0;
     else if ( !strcmp(key,"trace") ) {      else if ( !strcmp(key,"trace") ) {
            if ( value ) {        if ( value ) {
                u = BDY((LIST)value);          u = BDY((LIST)value);
            nd_nzlist = BDY((LIST)ARG2(u));          nd_nzlist = BDY((LIST)ARG2(u));
            nd_bpe = ZTOS((Q)ARG3(u));          nd_bpe = ZTOS((Q)ARG3(u));
            }        }
     } else if ( !strcmp(key,"f4red") ) {      } else if ( !strcmp(key,"f4red") ) {
        nd_f4red = ZTOS((Q)value);        nd_f4red = ZTOS((Q)value);
     } else if ( !strcmp(key,"rank0") ) {      } else if ( !strcmp(key,"rank0") ) {
             nd_rank0 = value?1:0;        nd_rank0 = value?1:0;
     } else if ( !strcmp(key,"splist") ) {      } else if ( !strcmp(key,"splist") ) {
             nd_splist = value?1:0;        nd_splist = value?1:0;
     } else if ( !strcmp(key,"check_splist") ) {      } else if ( !strcmp(key,"check_splist") ) {
       nd_check_splist = BDY((LIST)value);        nd_check_splist = BDY((LIST)value);
     } else if ( !strcmp(key,"sugarweight") ) {      } else if ( !strcmp(key,"sugarweight") ) {
       u = BDY((LIST)value);        u = BDY((LIST)value);
             n = length(u);        n = length(u);
             nd_sugarweight = MALLOC(n*sizeof(int));        nd_sugarweight = MALLOC(n*sizeof(int));
       for ( i = 0; i < n; i++, u = NEXT(u) )        for ( i = 0; i < n; i++, u = NEXT(u) )
                 nd_sugarweight[i] = ZTOS((Q)BDY(u));          nd_sugarweight[i] = ZTOS((Q)BDY(u));
       } else if ( !strcmp(key,"f4_td") ) {
         nd_f4_td = value?1:0;
       } else if ( !strcmp(key,"sba_f4step") ) {
         nd_sba_f4step = value?ZTOS((Q)value):0;
       } else if ( !strcmp(key,"sba_pot") ) {
         nd_sba_pot = value?1:0;
       } else if ( !strcmp(key,"sba_largelcm") ) {
         nd_sba_largelcm = value?1:0;
     }      }
     }    }
 }  }
   
 ND mdptond(DP d);  ND mdptond(DP d);
Line 9136  void nd_f4_lf_trace(LIST f,LIST v,int trace,int homo,s
Line 10165  void nd_f4_lf_trace(LIST f,LIST v,int trace,int homo,s
             ndv_homogenize((NDV)BDY(t),obpe,oadv,oepos,ompos);              ndv_homogenize((NDV)BDY(t),obpe,oadv,oepos,ompos);
     }      }
     if ( MaxDeg > 0 ) nocheck = 1;      if ( MaxDeg > 0 ) nocheck = 1;
     ret = ndv_setup(-2,m,fd0,nd_gbblock?1:0,0);      ret = ndv_setup(-2,m,fd0,nd_gbblock?1:0,0,0);
     if ( ret )      if ( ret )
       cand = nd_f4_lf_trace_main(m,&perm);        cand = nd_f4_lf_trace_main(m,&perm);
     if ( !ret || !cand ) {      if ( !ret || !cand ) {
Line 9252  NODE nd_f4_lf_trace_main(int m,int **indp)
Line 10281  NODE nd_f4_lf_trace_main(int m,int **indp)
             if ( DL_COMPARE(HDL(nfv),HDL(nfqv)) ) return 0;              if ( DL_COMPARE(HDL(nfv),HDL(nfqv)) ) return 0;
             ndv_removecont(m,nfv);              ndv_removecont(m,nfv);
             ndv_removecont(-2,nfqv);              ndv_removecont(-2,nfqv);
             nh = ndv_newps(-2,nfv,nfqv,1);              nh = ndv_newps(-2,nfv,nfqv);
             d = update_pairs(d,g,nh,0);              d = update_pairs(d,g,nh,0);
             g = update_base(g,nh);              g = update_base(g,nh);
         }          }
Line 9292  NDV vect64_to_ndv(mp_limb_t *vect,int spcol,int col,in
Line 10321  NDV vect64_to_ndv(mp_limb_t *vect,int spcol,int col,in
     }      }
 }  }
   
   NDV vect64_to_ndv_s(mp_limb_t *vect,int col,UINT *s0vect)
   {
       int j,k,len;
       UINT *p;
       UINT c;
       NDV r;
       NMV mr0,mr;
   
       for ( j = 0, len = 0; j < col; j++ ) if ( vect[j] ) len++;
       if ( !len ) return 0;
       else {
           mr0 = (NMV)MALLOC_ATOMIC_IGNORE_OFF_PAGE(nmv_adv*len);
           mr = mr0;
           p = s0vect;
           for ( j = k = 0; j < col; j++, p += nd_wpd )
             if ( (c = (UINT)vect[k++]) != 0 ) {
               ndl_copy(p,DL(mr)); CM(mr) = c; NMV_ADV(mr);
             }
           MKNDV(nd_nvar,mr0,len,r);
           return r;
       }
   }
   
 int nd_to_vect64(int mod,UINT *s0,int n,ND d,mp_limb_t *r)  int nd_to_vect64(int mod,UINT *s0,int n,ND d,mp_limb_t *r)
 {  {
     NM m;      NM m;
Line 9321  int nd_to_vect64(int mod,UINT *s0,int n,ND d,mp_limb_t
Line 10373  int nd_to_vect64(int mod,UINT *s0,int n,ND d,mp_limb_t
   
 #define MOD128(a,c,m) ((a)=(((c)!=0||((a)>=(m)))?(((((U128)(c))<<64)+(a))%(m)):(a)))  #define MOD128(a,c,m) ((a)=(((c)!=0||((a)>=(m)))?(((((U128)(c))<<64)+(a))%(m)):(a)))
   
 int ndv_reduce_vect64(int m,mp_limb_t *svect,mp_limb_t *cvect,int col,IndArray *imat,NM_ind_pair *rp0,int nred)  int ndv_reduce_vect64(int m,mp_limb_t *svect,mp_limb_t *cvect,int col,IndArray *imat,NM_ind_pair *rp0,int nred,SIG sig)
 {  {
     int i,j,k,len,pos,prev;      int i,j,k,len,pos,prev;
     mp_limb_t a,c,c1,c2;      mp_limb_t a,c,c1,c2;
Line 9342  int ndv_reduce_vect64(int m,mp_limb_t *svect,mp_limb_t
Line 10394  int ndv_reduce_vect64(int m,mp_limb_t *svect,mp_limb_t
         a = svect[k]; c = cvect[k];          a = svect[k]; c = cvect[k];
         MOD128(a,c,m);          MOD128(a,c,m);
         svect[k] = a; cvect[k] = 0;          svect[k] = a; cvect[k] = 0;
         if ( (c = svect[k]) != 0 ) {          if ( (c = svect[k]) != 0 && (sig == 0 || comp_sig(sig,rp0[i]->sig) > 0 ) ) {
             Nf4_red++;              Nf4_red++;
             maxrs = MAX(maxrs,rp0[i]->sugar);              maxrs = MAX(maxrs,rp0[i]->sugar);
             c = m-c; redv = nd_ps[rp0[i]->index];              c = m-c; redv = nd_ps[rp0[i]->index];
Line 9416  NODE nd_f4_red_mod64_main(int m,ND_pairs sp0,int nsp,U
Line 10468  NODE nd_f4_red_mod64_main(int m,ND_pairs sp0,int nsp,U
         nd_sp(m,0,sp,&spol);          nd_sp(m,0,sp,&spol);
         if ( !spol ) continue;          if ( !spol ) continue;
         nd_to_vect64(m,s0vect,col,spol,svect);          nd_to_vect64(m,s0vect,col,spol,svect);
         maxrs = ndv_reduce_vect64(m,svect,cvect,col,imat,rvect,nred);          maxrs = ndv_reduce_vect64(m,svect,cvect,col,imat,rvect,nred,0);
         for ( i = 0; i < col; i++ ) if ( svect[i] ) break;          for ( i = 0; i < col; i++ ) if ( svect[i] ) break;
         if ( i < col ) {          if ( i < col ) {
             spmat[sprow] = v = (mp_limb_t *)MALLOC_ATOMIC(spcol*sizeof(mp_limb_t));              spmat[sprow] = v = (mp_limb_t *)MALLOC_ATOMIC(spcol*sizeof(mp_limb_t));
Line 9548  int nd_gauss_elim_mod64(mp_limb_t **mat,int *sugar,ND_
Line 10600  int nd_gauss_elim_mod64(mp_limb_t **mat,int *sugar,ND_
 }  }
 #endif  #endif
   
   int nd_gauss_elim_mod64_s(mp_limb_t **mat,int *sugar,ND_pairs *spactive,int row,int col,int md,int *colstat,SIG *sig)
   {
     int i,j,k,l,rank,s,imin;
     mp_limb_t inv;
     mp_limb_t a;
     UINT c;
     mp_limb_t *t,*pivot,*pk;
     UINT *ck;
     UINT **cmat;
     UINT *ct;
     ND_pairs pair;
     SIG sg;
     int *used;
   
     used = (int *)MALLOC(row*sizeof(int));
     cmat = (UINT **)MALLOC(row*sizeof(UINT *));
     for ( i = 0; i < row; i++ ) {
       cmat[i] = MALLOC_ATOMIC(col*sizeof(UINT));
       bzero(cmat[i],col*sizeof(UINT));
     }
   
     for ( j = 0; j < col; j++ ) {
       for ( i = 0; i < row; i++ ) {
         a = mat[i][j]; c = cmat[i][j];
         MOD128(a,c,md);
         mat[i][j] = a; cmat[i][j] = 0;
       }
       for ( i = 0; i < row; i++ )
         if ( !used[i] && mat[i][j] ) break;
       if ( i == row ) {
         colstat[j] = 0;
         continue;
       } else {
         colstat[j] = 1;
         used[i] = 1;
       }
       /* column j is normalized */
       s = sugar[i];
       inv = invm((UINT)mat[i][j],md);
       /* normalize pivot row */
       for ( k = j, pk = mat[i]+j, ck = cmat[i]+j; k < col; k++, pk++, ck++ ) {
         a = *pk; c = *ck; MOD128(a,c,md); *pk = (a*inv)%md; *ck = 0;
       }
       for ( k = i+1; k < row; k++ ) {
         if ( (a = mat[k][j]) != 0 ) {
           sugar[k] = MAX(sugar[k],s);
           red_by_vect64(md,mat[k]+j,cmat[k]+j,mat[i]+j,(int)(md-a),col-j);
           Nf4_red++;
         }
       }
     }
     rank = 0;
     for ( i = 0; i < row; i++ ) {
       for ( j = 0; j < col; j++ )
         if ( mat[i][j] ) break;
       if ( j == col ) sugar[i] = -1;
       else rank++;
     }
     for ( i = 0; i < row; i++ ) GCFREE(cmat[i]);
     GCFREE(cmat);
     return rank;
   }
   
   NODE nd_f4_red_mod64_main_s(int m,ND_pairs sp0,int nsp,UINT *s0vect,int col,
           NM_ind_pair *rvect,int *rhead,IndArray *imat,int nred,NODE *syzlistp)
   {
       int spcol,sprow,a;
       int i,j,k,l,rank;
       NODE r0,r;
       ND_pairs sp;
       ND spol;
       mp_limb_t **spmat;
       mp_limb_t *svect,*cvect;
       mp_limb_t *v;
       int *colstat;
       struct oEGT eg0,eg1,eg2,eg_f4,eg_f4_1,eg_f4_2;
       int maxrs;
       int *spsugar;
       ND_pairs *spactive;
       SIG *spsig;
   
       get_eg(&eg0);
       /* elimination (1st step) */
       spmat = (mp_limb_t **)MALLOC(nsp*sizeof(mp_limb_t *));
       cvect = (mp_limb_t *)MALLOC(col*sizeof(mp_limb_t));
       spsugar = (int *)MALLOC(nsp*sizeof(int));
       spsig = (SIG *)MALLOC(nsp*sizeof(SIG));
       for ( a = sprow = 0, sp = sp0; a < nsp; a++, sp = NEXT(sp) ) {
           nd_sp(m,0,sp,&spol);
           if ( !spol ) {
             syzlistp[sp->sig->pos] = insert_sig(syzlistp[sp->sig->pos],sp->sig);
             continue;
           }
           svect = (mp_limb_t *)MALLOC(col*sizeof(mp_limb_t));
           nd_to_vect64(m,s0vect,col,spol,svect);
           maxrs = ndv_reduce_vect64(m,svect,cvect,col,imat,rvect,nred,spol->sig);
           for ( i = 0; i < col; i++ ) if ( svect[i] ) break;
           if ( i < col ) {
               spmat[sprow] = svect;
               spsugar[sprow] = MAX(maxrs,SG(spol));
               spsig[sprow] = sp->sig;
               sprow++;
           } else {
             syzlistp[sp->sig->pos] = insert_sig(syzlistp[sp->sig->pos],sp->sig);
           }
           nd_free(spol);
       }
       get_eg(&eg1); init_eg(&eg_f4_1); add_eg(&eg_f4_1,&eg0,&eg1); add_eg(&f4_elim1,&eg0,&eg1);
       if ( DP_Print ) {
           fprintf(asir_out,"elim1=%.3fsec,",eg_f4_1.exectime);
           fflush(asir_out);
       }
       /* free index arrays */
       for ( i = 0; i < nred; i++ ) GCFREE(imat[i]->index.c);
   
       /* elimination (2nd step) */
       colstat = (int *)MALLOC(col*sizeof(int));
       rank = nd_gauss_elim_mod64_s(spmat,spsugar,0,sprow,col,m,colstat,spsig);
       r0 = 0;
       for ( i = 0; i < sprow; i++ ) {
           if ( spsugar[i] >= 0 ) {
             NEXTNODE(r0,r);
             BDY(r) = vect64_to_ndv_s(spmat[i],col,s0vect);
             SG((NDV)BDY(r)) = spsugar[i];
             ((NDV)BDY(r))->sig = spsig[i];
           } else
             syzlistp[spsig[i]->pos] = insert_sig(syzlistp[spsig[i]->pos],spsig[i]);
           GCFREE(spmat[i]);
       }
       if ( r0 ) NEXT(r) = 0;
       get_eg(&eg2); init_eg(&eg_f4_2); add_eg(&eg_f4_2,&eg1,&eg2); add_eg(&f4_elim2,&eg1,&eg2);
       init_eg(&eg_f4); add_eg(&eg_f4,&eg0,&eg2);
       if ( DP_Print ) {
           fprintf(asir_out,"elim2=%.3fsec,",eg_f4_2.exectime);
           fprintf(asir_out,"nsp=%d,nred=%d,spmat=(%d,%d),rank=%d ",
               nsp,nred,sprow,col,rank);
           fprintf(asir_out,"%.3fsec,",eg_f4.exectime);
       }
       return r0;
   }
   
   NODE nd_f4_red_s(int m,ND_pairs sp0,int trace,UINT *s0vect,int col,NODE rp0,NODE *syzlistp)
   {
     IndArray *imat;
     int nsp,nred,i,start;
     int *rhead;
     NODE r0,rp;
     ND_pairs sp;
     NM_ind_pair *rvect;
     UINT *s;
     int *s0hash;
     struct oEGT eg0,eg1,eg_conv;
   
     for ( sp = sp0, nsp = 0; sp; sp = NEXT(sp), nsp++ );
     nred = length(rp0);
     imat = (IndArray *)MALLOC(nred*sizeof(IndArray));
     rhead = (int *)MALLOC(col*sizeof(int));
     for ( i = 0; i < col; i++ ) rhead[i] = 0;
   
     /* construction of index arrays */
     get_eg(&eg0);
     if ( DP_Print ) {
       fprintf(asir_out,"%dx%d,",nsp+nred,col);
       fflush(asir_out);
     }
     rvect = (NM_ind_pair *)MALLOC(nred*sizeof(NM_ind_pair));
     for ( start = 0, rp = rp0, i = 0; rp; i++, rp = NEXT(rp) ) {
       rvect[i] = (NM_ind_pair)BDY(rp);
       imat[i] = nm_ind_pair_to_vect_compress(trace,s0vect,col,rvect[i],start);
       rhead[imat[i]->head] = 1;
       start = imat[i]->head;
     }
     get_eg(&eg1); init_eg(&eg_conv); add_eg(&eg_conv,&eg0,&eg1); add_eg(&f4_conv,&eg0,&eg1);
     if ( DP_Print ) {
       fprintf(asir_out,"conv=%.3fsec,",eg_conv.exectime);
       fflush(asir_out);
     }
     if ( m > 0 )
       r0 = nd_f4_red_mod64_main_s(m,sp0,nsp,s0vect,col,rvect,rhead,imat,nred,syzlistp);
     else
   //    r0 = nd_f4_red_q_main_s(sp0,nsp,trace,s0vect,col,rvect,rhead,imat,nred);
       error("nd_f4_red_q_main_s : not implemented yet");
     return r0;
   }
   
   INLINE int ndl_find_reducer_minsig(UINT *dg)
   {
     RHist r;
     int i,singular,ret,d,k,imin;
     SIG t;
     static int wpd,nvar;
     static SIG quo,quomin;
     static UINT *tmp;
   
     if ( !quo || nvar != nd_nvar ) { NEWSIG(quo); NEWSIG(quomin); }
     if ( wpd != nd_wpd ) {
       wpd = nd_wpd;
       tmp = (UINT *)MALLOC(wpd*sizeof(UINT));
     }
   #if 0
     d = ndl_hash_value(dg);
     for ( r = nd_red[d], k = 0; r; r = NEXT(r), k++ ) {
       if ( ndl_equal(dg,DL(r)) ) {
         return r->index;
       }
     }
   #endif
     imin = -1;
     for ( i = 0; i < nd_psn; i++ ) {
       r = nd_psh[i];
       if ( ndl_reducible(dg,DL(r)) ) {
         ndl_sub(dg,DL(r),tmp);
         _ndltodl(tmp,DL(quo));
         _addtodl(nd_nvar,DL(nd_psh[i]->sig),DL(quo));
         quo->pos = nd_psh[i]->sig->pos;
         if ( imin < 0 || comp_sig(quomin,quo) > 0 ) {
           t = quo; quo = quomin; quomin = t;
           imin = i;
         }
       }
     }
     if ( imin == -1 ) return nd_psn;
     else {
   #if 0
       nd_append_red(dg,i);
   #endif
       return imin;
     }
   }
   
   int nd_symbolic_preproc_s(PGeoBucket bucket,int trace,UINT **s0vect,NODE *r)
   {
     NODE rp0,rp;
     NM mul,head,s0,s;
     int index,col,i,sugar;
     RHist h;
     UINT *s0v,*p;
     NM_ind_pair pair;
     ND red;
     NDV *ps;
     SIG sig;
   
     s0 = 0; rp0 = 0; col = 0;
     if ( nd_demand )
       ps = trace?nd_ps_trace_sym:nd_ps_sym;
     else
       ps = trace?nd_ps_trace:nd_ps;
     while ( 1 ) {
       head = remove_head_pbucket_symbolic(bucket);
       if ( !head ) break;
       if ( !s0 ) s0 = head;
       else NEXT(s) = head;
       s = head;
       index = ndl_find_reducer_minsig(DL(head));
       if ( index >= 0 && index < nd_psn ) {
         h = nd_psh[index];
         NEWNM(mul);
         ndl_sub(DL(head),DL(h),DL(mul));
         if ( ndl_check_bound2(index,DL(mul)) )
           return 0;
         sugar = TD(DL(mul))+SG(ps[index]);
         NEWSIG(sig);
         _ndltodl(DL(mul),DL(sig));
         _addtodl(nd_nvar,DL(nd_psh[index]->sig),DL(sig));
         sig->pos = nd_psh[index]->sig->pos;
         MKNM_ind_pair(pair,mul,index,sugar,sig);
         red = ndv_mul_nm_symbolic(mul,ps[index]);
         add_pbucket_symbolic(bucket,nd_remove_head(red));
         NEXTNODE(rp0,rp); BDY(rp) = (pointer)pair;
       }
       col++;
     }
     if ( rp0 ) NEXT(rp) = 0;
     NEXT(s) = 0;
     s0v = (UINT *)MALLOC_ATOMIC(col*nd_wpd*sizeof(UINT));
     for ( i = 0, p = s0v, s = s0; i < col;
       i++, p += nd_wpd, s = NEXT(s) ) ndl_copy(DL(s),p);
     *s0vect = s0v;
     *r = rp0;
   
     return col;
   }
   
   NODE nd_sba_f4(int m,int **indp)
   {
     int i,nh,stat,index,f4red,f4step;
     int col,rank,len,k,j,a,sugar,nbase,psugar,ms;
     NODE r,g,rp0,nflist;
     ND_pairs d,l,t;
     ND h,nf;
     NDV nfv;
     union oNDC hc;
     UINT *s0vect;
     UINT c;
     PGeoBucket bucket;
     NODE *syzlist;
     SIG sig;
     struct oEGT eg0,eg1,eg_f4;
     struct oEGT eg2,eg_update,eg_remove,eg_large,eg_nf,eg_nfzero;
   
     Nf4_red=0;
     d = 0;
     syzlist = (NODE *)MALLOC(nd_psn*sizeof(NODE));
     for ( i = 0; i < nd_psn; i++ ) {
       d = update_pairs_s(d,i,syzlist);
     }
     nd_nbase = nd_psn;
     f4red = 1;
     psugar = 0;
     f4step = 0;
     while ( d ) {
       for ( t = d, ms = SG(d); t; t = NEXT(t) )
         if ( SG(t) < ms ) ms = SG(t);
       if ( ms == psugar && f4step >= nd_sba_f4step ) {
   again:
         l = d; d = d->next;
         if ( small_lcm(l) ) {
           if ( DP_Print ) fprintf(asir_out,"M");
           continue;
         }
         sig = l->sig;
         stat = nd_sp(m,0,l,&h);
         if ( !stat ) {
           NEXT(l) = d; d = l;
           d = nd_reconstruct(0,d);
           goto again;
         }
     get_eg(&eg1);
     #if USE_GEOBUCKET
         stat = m?nd_nf_pbucket_s(m,h,nd_ps,!Top,&nf):nd_nf_s(m,0,h,nd_ps,!Top,&nf);
     #else
         stat = nd_nf_s(m,0,h,nd_ps,!Top,&nf);
     #endif
     get_eg(&eg2);
         if ( !stat ) {
           NEXT(l) = d; d = l;
           d = nd_reconstruct(0,d);
           goto again;
         } else if ( stat == -1 ) {
           if ( DP_Print ) { printf("S"); fflush(stdout); }
           FREENDP(l);
         } else if ( nf ) {
           if ( DP_Print ) { printf("+"); fflush(stdout); }
           add_eg(&eg_nf,&eg1,&eg2);
           hc = HCU(nf);
           nd_removecont(m,nf);
           nfv = ndtondv(m,nf); nd_free(nf);
           nh = ndv_newps(m,nfv,0);
   
           d = update_pairs_s(d,nh,syzlist);
           nd_sba_pos[sig->pos] = append_one(nd_sba_pos[sig->pos],nh);
           FREENDP(l);
         } else {
           add_eg(&eg_nfzero,&eg1,&eg2);
          // syzygy
     get_eg(&eg1);
           d = remove_spair_s(d,sig);
     get_eg(&eg2); add_eg(&eg_remove,&eg1,&eg2);
           syzlist[sig->pos] = insert_sig(syzlist[sig->pos],sig);
           if ( DP_Print ) { printf("."); fflush(stdout); }
           FREENDP(l);
         }
       } else {
         if ( ms != psugar ) f4step = 1;
         else f4step++;
   again2:
         psugar = ms;
         l = nd_minsugarp_s(d,&d);
         sugar = nd_sugarweight?d->sugar2:SG(d);
         bucket = create_pbucket();
         stat = nd_sp_f4(m,0,l,bucket);
         if ( !stat ) {
           for ( t = l; NEXT(t); t = NEXT(t) );
           NEXT(t) = d; d = l;
           d = nd_reconstruct(0,d);
           goto again2;
         }
         if ( bucket->m < 0 ) continue;
         col = nd_symbolic_preproc_s(bucket,0,&s0vect,&rp0);
         if ( !col ) {
           for ( t = l; NEXT(t); t = NEXT(t) )
             ;
           NEXT(t) = d; d = l;
           d = nd_reconstruct(0,d);
           goto again2;
         }
         if ( DP_Print ) fprintf(asir_out,"\nsugar=%d,",psugar);
         nflist = nd_f4_red_s(m,l,0,s0vect,col,rp0,syzlist);
         /* adding new bases */
         for ( r = nflist; r; r = NEXT(r) ) {
           nfv = (NDV)BDY(r);
           if ( nd_f4_td ) SG(nfv) = nd_tdeg(nfv);
           ndv_removecont(m,nfv);
           nh = ndv_newps(m,nfv,0);
           d = update_pairs_s(d,nh,syzlist);
           nd_sba_pos[nfv->sig->pos] = append_one(nd_sba_pos[nfv->sig->pos],nh);
         }
         for ( i = 0; i < nd_nbase; i++ )
           for ( r = syzlist[i]; r; r = NEXT(r) )
               d = remove_spair_s(d,(SIG)BDY(r));
         d = remove_large_lcm(d);
         if ( DP_Print ) {
           fprintf(asir_out,"f4red=%d,gblen=%d",f4red,nd_psn); fflush(asir_out);
         }
         f4red++;
       }
     }
     if ( DP_Print ) {
       fprintf(asir_out,"\nnumber of red=%d,",Nf4_red);
     }
     g = conv_ilist_s(nd_demand,0,indp);
     return g;
   }
Legend:
Removed from v.1.21  
changed lines
  Added in v.1.38
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>