[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.55 and 1.59

version 1.55, 2021/11/29 09:19:33 version 1.59, 2021/12/14 01:02:40
Line 1 
Line 1 
 /* $OpenXM: OpenXM_contrib2/asir2018/engine/nd.c,v 1.54 2021/10/29 20:50:02 noro Exp $ */  /* $OpenXM: OpenXM_contrib2/asir2018/engine/nd.c,v 1.58 2021/12/07 08:58:08 noro Exp $ */
   
 #include "nd.h"  #include "nd.h"
   
Line 2479  LIST compute_splist()
Line 2479  LIST compute_splist()
 }  }
   
 typedef struct oHPDATA {  typedef struct oHPDATA {
   P hn; // HP(t)=hn(t)/(1-t)^n    int n;
   int len;    P hn; // HP(t)=hn(t)/((1-t^w0)*...*(1-t^w(n-1)))
   P *head; // hp(i)=head[i] (i=0,...,len-1)  
   P hp; // dim Hm(i)=hp(i) (i >= len)  
   VECT x; // BDY(x)[i] = <<0,...,1,...,0>>    VECT x; // BDY(x)[i] = <<0,...,1,...,0>>
   P *plist; // plist[i]=(1-t)^i    int *w;
 } *HPDATA;  } *HPDATA;
   
 void make_reduced(VECT b,int nv);  void make_reduced(VECT b,int nv);
 void mhp_rec(VECT b,VECT x,P t,P *r);  void mhp_rec(VECT b,VECT x,P t,P *r);
 P mhp_ctop(P *r,P *plist,int n);  P mhp_ctop(P *r,P *plist,int n);
 void mhp_to_hf(VL vl,P hp,int n,P *plist,VECT *head,P *hf);  
 DL monomial_colon(DL a,DL b,int n);  DL monomial_colon(DL a,DL b,int n);
 LIST dp_monomial_hilbert_poincare(VECT b,VECT x,P *plist);  LIST dp_monomial_hilbert_poincare(VECT b,VECT x);
   LIST dp_monomial_hilbert_poincare_weight(VECT b,VECT x,int *w);
   
 int hpvalue(HPDATA data,int d)  
 {  
   P *head;  
   int len;  
   P hp,val;  
   Z dz;  
   
   head = data->head;  
   len = data->len;  
   hp = data->hp;  
   if ( d < len )  
     return ZTOS((Z)head[d]);  
   else {  
     STOZ(d,dz);  
     substp(CO,hp,hp->v,(P)dz,&val);  
     return ZTOS((Z)val);  
   }  
 }  
   
 void setup_hpdata(HPDATA final,HPDATA current)  void setup_hpdata(HPDATA final,HPDATA current)
 {  {
   int n,i;    int n,i,wlen;
   P *r;  
   DL *p;    DL *p;
   P tv;    VECT b,x;
   VECT b,x,head;  
   DL dl;    DL dl;
     LIST weight;
     LIST ret;
     int *w;
     NODE nd;
   
   n = nd_nvar;    final->n = n = nd_nvar;
   final->hn = (P)ARG0(nd_hpdata);    final->hn = (P)BDY(nd_hpdata);
   head = (VECT)ARG2(nd_hpdata);  #if 0
   final->len = head->len;    if ( NEXT(nd_hpdata) != 0 && (weight=(LIST)BDY(NEXT(nd_hpdata))) != 0 ) {
   final->head = (P *)BDY(head);      wlen = length(BDY(weight));
   final->hp = (P)ARG3(nd_hpdata);      if ( n != wlen )
   final->plist = (P *)BDY((VECT)ARG4(nd_hpdata));        error("setup_hpdata : 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
       w = 0;
   #else
     w = current_dl_weight_vector;
   #endif
   MKVECT(x,n);    MKVECT(x,n);
   for ( i = 0; i < n; i++ ) {    for ( i = 0; i < n; i++ ) {
     NEWDL(dl,n); dl->d[i] = 1; dl->td = 1; BDY(x)[i] = dl;      NEWDL(dl,n); dl->d[i] = 1; dl->td = 1; BDY(x)[i] = dl;
   }    }
   final->x = x;    final->x = x;
     final->w = w;
   
   r = (P *)CALLOC(n+1,sizeof(P));  
   MKVECT(b,nd_psn); p = (DL *)BDY(b);    MKVECT(b,nd_psn); p = (DL *)BDY(b);
   for ( i = 0; i < nd_psn; i++ ) {    for ( i = 0; i < nd_psn; i++ ) {
     p[i] = ndltodl(n,nd_psh[i]->dl);      p[i] = ndltodl(n,nd_psh[i]->dl);
   }    }
   make_reduced(b,n);    if ( w ) {
   makevar("t",&tv);      ret = dp_monomial_hilbert_poincare_weight(b,x,w);
   mhp_rec(b,x,tv,r);    } else
   current->hn = mhp_ctop(r,final->plist,n);      ret = dp_monomial_hilbert_poincare(b,x);
   mhp_to_hf(CO,current->hn,n,final->plist,&head,&current->hp);    current->n = n;
   current->head = (P *)BDY(head);    current->hn = (P)BDY(BDY(ret));
   current->len = head->len;  
   current->x = x;    current->x = x;
   current->plist = final->plist;    current->w = w;
 }  }
   
 void update_hpdata(HPDATA current,int nh,int do_hf)  int comp_hn(P a, P b)
 {  {
     P s;
     DCP dc;
   
     subp(CO,a,b,&s);
     if ( !s ) return -1;
     else if ( OID(s) == 1 ) return 0;
     else {
       for ( dc = DC(s); NEXT(dc); dc = NEXT(dc) );
       return (int)ZTOS((Z)dc->d);
     }
   }
   
   void update_hpdata(HPDATA current,int nh)
   {
   NODE data1,nd,t;    NODE data1,nd,t;
   DL new,dl;    DL new,dl;
   int len,i,n;    int len,i,n;
Line 2569  void update_hpdata(HPDATA current,int nh,int do_hf)
Line 2571  void update_hpdata(HPDATA current,int nh,int do_hf)
     p[i] = monomial_colon(ndltodl(n,nd_psh[i]->dl),new,n);      p[i] = monomial_colon(ndltodl(n,nd_psh[i]->dl),new,n);
   }    }
   // compute HP(I:new)    // compute HP(I:new)
   list1 = dp_monomial_hilbert_poincare(b,current->x,current->plist);    if ( current->w )
       list1 = dp_monomial_hilbert_poincare_weight(b,current->x,current->w);
     else
       list1 = dp_monomial_hilbert_poincare(b,current->x);
   data1 = BDY((LIST)list1);    data1 = BDY((LIST)list1);
   // HP(I+<new>) = H(I)-t^d*H(I:new), d=tdeg(new)    // HP(I+<new>) = H(I)-t^d*H(I:new), d=tdeg(new)
   makevar("t",&tv); UTOZ(new->td,dz);    makevar("t",&tv); UTOZ(new->td,dz);
Line 2577  void update_hpdata(HPDATA current,int nh,int do_hf)
Line 2582  void update_hpdata(HPDATA current,int nh,int do_hf)
   mulp(CO,(P)ARG0(data1),td,&s);    mulp(CO,(P)ARG0(data1),td,&s);
   subp(CO,current->hn,s,&hn);    subp(CO,current->hn,s,&hn);
   current->hn = hn;    current->hn = hn;
   if ( do_hf ) {  
     mhp_to_hf(CO,hn,n,current->plist,&head,&hpoly);  
     current->head = (P *)BDY(head);  
     current->len = head->len;  
     current->hp = hpoly;  
   }  
 }  }
   
 ND_pairs nd_remove_same_sugar( ND_pairs d, int sugar)  ND_pairs nd_remove_same_sugar( ND_pairs d, int sugar)
Line 2651  again:
Line 2650  again:
       }        }
       sugar = SG(l);        sugar = SG(l);
       if ( DP_Print ) fprintf(asir_out,"%d",sugar);        if ( DP_Print ) fprintf(asir_out,"%d",sugar);
       if ( nd_hpdata ) {  
         if ( !compp(CO,final_hpdata.hn,current_hpdata.hn) )  
           break;  
         else {  
           final_hpvalue = hpvalue(&final_hpdata,sugar);  
           if ( final_hpvalue == hpvalue(&current_hpdata,sugar) ) {  
 //            if ( DP_Print ) fprintf(asir_out,"done.\n",sugar);  
             d = nd_remove_same_sugar(d,sugar);  
             continue;  
           }  
         }  
       }  
     }      }
     stat = nd_sp(m,0,l,&h);      stat = nd_sp(m,0,l,&h);
     if ( !stat ) {      if ( !stat ) {
Line 2721  get_eg(&eg2); add_eg(&eg_update,&eg1,&eg2);
Line 2708  get_eg(&eg2); add_eg(&eg_update,&eg1,&eg2);
       g = update_base(g,nh);        g = update_base(g,nh);
       FREENDP(l);        FREENDP(l);
       if ( nd_hpdata ) {        if ( nd_hpdata ) {
         update_hpdata(&current_hpdata,nh,1);          int dg,sugar0;
         if ( final_hpvalue == hpvalue(&current_hpdata,sugar) ) {  
 //          if ( DP_Print ) fprintf(asir_out,"sugar=%d done.\n",sugar);          update_hpdata(&current_hpdata,nh);
           d = nd_remove_same_sugar(d,sugar);          dg = comp_hn(final_hpdata.hn,current_hpdata.hn);
           if ( dg < 0 ) {
              fprintf(asir_out,"We found a gb\n");
              d = 0;
         }          }
           sugar0 = sugar;
           while ( d && dg > sugar0 ) {
             if ( DP_Print ) { fprintf(asir_out,"[%d]",sugar0); fflush(asir_out); }
             d = nd_remove_same_sugar(d,sugar0);
             sugar0++;
           }
       }        }
     } else {      } else {
       Nnfz++;        Nnfz++;
Line 2736  get_eg(&eg2); add_eg(&eg_update,&eg1,&eg2);
Line 2732  get_eg(&eg2); add_eg(&eg_update,&eg1,&eg2);
         MKNODE(t,list,nd_alltracelist);          MKNODE(t,list,nd_alltracelist);
         nd_alltracelist = t; nd_tracelist = 0;          nd_alltracelist = t; nd_tracelist = 0;
       }        }
       if ( DP_Print ) { printf("."); fflush(stdout); }        if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); }
         FREENDP(l);          FREENDP(l);
     }      }
   }    }
   conv_ilist(nd_demand,0,g,indp);    conv_ilist(nd_demand,0,g,indp);
   if ( !checkonly && DP_Print ) {    if ( !checkonly && DP_Print ) {
     printf("\nnd_gb done. Nnd_add=%d,Npairs=%d, Nnfnz=%d,Nnfz=%d,",Nnd_add,Npairs,Nnfnz,Nnfz);      fprintf(asir_out,"\nnd_gb done. Nnd_add=%d,Npairs=%d, Nnfnz=%d,Nnfz=%d,",Nnd_add,Npairs,Nnfnz,Nnfz);
     printf("Nremoved=%d\n",NcriB+NcriMF+Ncri2);      fprintf(asir_out,"Nremoved=%d\n",NcriB+NcriMF+Ncri2);
     fflush(asir_out);      fflush(asir_out);
   }    }
   if ( DP_Print ) {    if ( DP_Print ) {
Line 3064  init_eg(&eg_remove);
Line 3060  init_eg(&eg_remove);
   Nnominimal = 0;    Nnominimal = 0;
   Nredundant = 0;    Nredundant = 0;
   ngen = nd_psn;    ngen = nd_psn;
   if ( !do_weyl ) {    if ( !do_weyl || nd_sba_inputisgb ) {
     for ( i = 0; i < nd_psn; i++ )      for ( i = 0; i < nd_psn; i++ )
       for ( j = i+1; j < nd_psn; j++ ) {        for ( j = i+1; j < nd_psn; j++ ) {
         sig = trivial_sig(i,j);          sig = trivial_sig(i,j);
Line 3172  get_eg(&eg2); 
Line 3168  get_eg(&eg2); 
        nd_sba_pos[sig->pos] = append_one(nd_sba_pos[sig->pos],nh);         nd_sba_pos[sig->pos] = append_one(nd_sba_pos[sig->pos],nh);
        if ( nd_hpdata ) {         if ( nd_hpdata ) {
          get_eg(&eg1);           get_eg(&eg1);
          update_hpdata(&current_hpdata,nh,0);           update_hpdata(&current_hpdata,nh);
          get_eg(&eg2); add_eg(&eg_hpdata,&eg1,&eg2);           get_eg(&eg2); add_eg(&eg_hpdata,&eg1,&eg2);
          if ( !compp(CO,final_hpdata.hn,current_hpdata.hn) ) {           if ( !compp(CO,final_hpdata.hn,current_hpdata.hn) ) {
            if ( DP_Print ) { printf("\nWe found a gb.\n"); }             if ( DP_Print ) { printf("\nWe found a gb.\n"); }
Line 3441  again:
Line 3437  again:
 #endif  #endif
       sugar = SG(l);        sugar = SG(l);
       if ( DP_Print ) fprintf(asir_out,"%d",sugar);        if ( DP_Print ) fprintf(asir_out,"%d",sugar);
       if ( nd_hpdata ) {  
         if ( !compp(CO,final_hpdata.hn,current_hpdata.hn) )  
           break;  
         else {  
           final_hpvalue = hpvalue(&final_hpdata,sugar);  
           if ( final_hpvalue == hpvalue(&current_hpdata,sugar) ) {  
 //            if ( DP_Print ) fprintf(asir_out,"sugar=%d done.\n",sugar);  
             d = nd_remove_same_sugar(d,sugar);  
             continue;  
           }  
         }  
       }  
     }      }
     stat = nd_sp(m,0,l,&h);      stat = nd_sp(m,0,l,&h);
     if ( !stat ) {      if ( !stat ) {
Line 3527  again:
Line 3511  again:
         d = update_pairs(d,g,nh,0);          d = update_pairs(d,g,nh,0);
         g = update_base(g,nh);          g = update_base(g,nh);
         if ( nd_hpdata ) {          if ( nd_hpdata ) {
           update_hpdata(&current_hpdata,nh,1);            int dg,sugar0;
           if ( final_hpvalue == hpvalue(&current_hpdata,sugar) ) {  
 //            if ( DP_Print ) fprintf(asir_out,"sugar=%d done.\n",sugar);            update_hpdata(&current_hpdata,nh);
             d = nd_remove_same_sugar(d,sugar);            dg = comp_hn(final_hpdata.hn,current_hpdata.hn);
             if ( dg < 0 ) {
                fprintf(asir_out,"We found a gb\n");
                d = 0;
           }            }
             sugar0 = sugar;
             while ( d && dg > sugar0 ) {
               if ( DP_Print ) { fprintf(asir_out,"[%d]",sugar0); fflush(asir_out); }
               d = nd_remove_same_sugar(d,sugar0);
               sugar0++;
             }
         }          }
       } else {        } else {
         if ( DP_Print ) { printf("*"); fflush(stdout); }          if ( DP_Print ) { fprintf(asir_out,"*"); fflush(asir_out); }
       }        }
     } else {      } else {
       if ( DP_Print ) { printf("."); fflush(stdout); }        if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); }
     }      }
     FREENDP(l);      FREENDP(l);
   }    }
Line 3549  again:
Line 3542  again:
     }      }
   }    }
   conv_ilist(nd_demand,1,g,indp);    conv_ilist(nd_demand,1,g,indp);
   if ( DP_Print ) { printf("\nnd_gb_trace done.\n"); fflush(stdout); }    if ( DP_Print ) { fprintf(asir_out,"\nnd_gb_trace done.\n"); fflush(asir_out); }
   return g;    return g;
 }  }
   
Legend:
Removed from v.1.55  
changed lines
  Added in v.1.59
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>