Logo Search packages:      
Sourcecode: palp version File versions  Download package

cws.c

#include "Global.h"
#include "LG.h"
#include "Rat.h"

FILE *inFILE, *outFILE;

#define     Only_IP_CWS                   1
#define TRANS_INFO_FOR_IP_WEIGHTS   0

#define NFmax  10             /* maximal number of WS-files */

#define SIMPLEX_POINT_Nmax          50

typedef struct {  int  u[NFmax];  int nu; }     CWS_type;

void  PrintCWSUsage(/* char *c */)
{     printf(
"%s%s%s%s%s%s\n%s\n\n%s\n%s\n%s\n%s\n%s\n%s\n%s\n%s\n\n%s\n%s\n%s\n%s\n%s\n",
"This is `cws.x': create weight systems and combined weight systems.\n",
"Usage:   cws.x -<options>;     the first option must be `w', `c', `i', or `h'"
,".\n\n",
"Options: -h        print this information\n\n",

"         -w# [L H] make IP weight systems for #-dimensional polytopes.\n",
"                   For #>4 the lowest and highest degrees L<=H are required.",
"             -r/-t make reflexive/transversal weight systems (optional).",

"         -c#       make combined weight systems for #-dimensional polytopes."
,"                   For #<=4 all relevant combinations are made by default,",
"                   otherwise the following option is required:",
"             -n[#] followed by the names wf_1 ... wf_# of weight files",
"                   currently #=2,3 are implemented.",
"              [-t] followed by # numbers n_i specifies the CWS-type, i.e.",
"                   the numbers n_i of weights to be selected from wf_i.",
"                   Currently all cases with n_i<=2 are implemented.",

"          -i       compute the polytope data M:p v [F:f] N:p [v] for all IP",
"                   CWS, where p and v denote the numbers of lattice points",
"                   and vertices of a dual pair of IP polytopes; an entry  ",
"                   F:f and no v for N indicates a non-reflexive `dual pair'.",
"              -f   use as filter; otherwise parameters denote I/O files\n");
      exit(0);
}

void Die(char *comment){ printf("\n%s\n",comment); exit(0);}

int Read_Weight(Weight *);

int READ_Weight(Weight *_W, FILE *INFILE){
  inFILE=INFILE;
  return Read_Weight(_W);
}

int READ_CWS_PP(CWS *_CW, PolyPointList *_P, FILE *INFILE){
  inFILE=INFILE;
  return Read_CWS_PP(_CW, _P);
}

void Print_CWS(CWS *_W);
void Conv(int narg, char* fn[]);
void SimplexPointCount(int narg, char* fn[]);
void Init_IP_Weights(int narg, char* fn[]);
void Init_IP_CWS(int narg, char* fn[]);
void IP_Poly_Data(int narg, char* fn[]);
void Make_CWS_Points(CWS *, PolyPointList *);
void Npoly2cws(int narg, char* fn[]);
void PrintCWSextUsage(char *c)
{    printf( "This is `%s': -x gives undocumented extensions:\n",c);
            puts("              -ip    printf PolyPointList");
            puts("              -id    printf dual PolyPointList");
            puts("              -N     make CWS for PPL in N lattice");
            puts("              -p#    [infile1] [infile2] makes cartesian product"
     );     puts("                     of Vertices. # dimensions are identified.");
            puts("              -S     count simplex points for weight system");
            puts("              -L     count using LattE (-> count redcheck cdd)");
}

int  main (int narg, char* fn[])
{    inFILE=stdin; outFILE=stdout;
     if(narg==1) {printf("\nFor help type `%s -h'\n\n",fn[0]);exit(0);}
     if((fn[1][0]!='-')||(fn[1][1]=='h')) PrintCWSUsage(/*fn[0]*/);
     else if(fn[1][1]=='w') Init_IP_Weights(narg, fn);
     else if(fn[1][1]=='c') Init_IP_CWS(narg, fn);
     else if(fn[1][1]=='i') IP_Poly_Data(narg, fn);
     else if(fn[1][1]=='N') Npoly2cws(narg, fn);
     else if(fn[1][1]=='p') Conv(narg, fn);
     else if(fn[1][1]=='x') PrintCWSextUsage(fn[0]);
     else if(fn[1][1]=='S') SimplexPointCount(narg, fn);
     else if(fn[1][1]=='L') SimplexPointCount(narg, fn);
     else PrintCWSUsage(/*fn[0]*/); return 0;
}
void Make_IP_Weights(int d, int Dmin, int Dmax,int rFlag,int tFlag);
void Make_34_Weights(int d, int tFlag);
void Init_IP_Weights(int narg, char* fn[])
{    int n=1,d,L=0,H=0,rf=0,tf=0; char *c=&fn[1][2];
     if(narg>2) if(c[0]==0) c=fn[++n];
     if(!IsDigit(c[0])){ puts("-w must be followed by a number");exit(0);}
     if(POLY_Dmax<(d=atoi(c))){printf("Increase POLY_Dmax to %d\n",d);exit(0);}
     if(++n<narg) if((fn[n][0]!='-')&&(IsDigit(fn[n][0])))
     {      L=atoi(fn[n]); assert(++n<narg); assert(IsDigit(fn[n][0]));
      H=atoi(fn[n]); assert(L<=H); n++;
     }      n--;
     while(++n<narg) if(fn[n][0]=='-')
     {  if(fn[n][1]=='r')rf=1; else if(fn[n][1]=='t')tf=1; 
      else {printf("Illegal option %s\n",fn[n]);exit(0);}
     }
     if(++n<narg) {printf("Want %s as output file?\n",fn[n]);exit(0);}
     if(n<narg) {puts("Too many arguments");exit(0);}
     if(H) Make_IP_Weights(d,L,H,rf,tf); else Make_34_Weights(d,tf);
     {  ;
     }
}

int  VP_2_CWS(Long *V[], int d, int v, CWS *W);
void Npoly2cws(int narg, char* fn[])
{    int n=2; CWS W; EqList E; VertexNumList V; Long *X[VERT_Nmax]; FILE *OF;
     PolyPointList *P=(PolyPointList *)malloc(sizeof(PolyPointList)); 
     assert(P!=NULL); assert(!strcmp(fn[1],"-N")); inFILE=stdin;outFILE=stdout;
     if(narg>2) {if(fn[2][0]=='-'){assert(fn[2][1]=='f');inFILE=NULL;} else
     {      inFILE=fopen(fn[2],"r"); assert(NULL!=inFILE);
        if(narg>3) {outFILE=fopen(fn[3],"w"); assert(NULL!=outFILE);}
     }}     OF=outFILE; while(Read_CWS_PP(&W,P))
     {  if(W.N) Die("Only PPL-input in Npoly2cws!");
      if(!IP_Check(P,&V,&E)) Die("Not IP!"); Sort_VL(&V);
      for(n=0;n<V.nv;n++) X[n]=P->x[V.v[n]];
      if(VP_2_CWS(X,P->n,V.nv,&W)) {Print_CWS(&W);fprintf(outFILE,"\n");}
      else {outFILE=stderr; Print_PPL(P,"CWS not found"); outFILE=OF;}
     }
}

void Make_IP_CWS(int narg, char* fn[]);
void Make_34_CWS(int d);
void Init_IP_CWS(int narg, char* fn[])
{    int d,n=1,nop=0; char *c=&fn[1][2]; if(narg>2) if(c[0]==0) c=fn[++n];
     if(!IsDigit(c[0])){ puts("-c must be followed by a number");exit(0);}
     if(POLY_Dmax<(d=atoi(c))){printf("Increase POLY_Dmax to %d\n",d);exit(0);}
     if(++n<narg) if((fn[n][0]=='-')&&(fn[n][1]=='n')) nop=1;
     if(nop) Make_IP_CWS(narg,fn); else if(d<=4) Make_34_CWS(d); else Die(
       "`-c#' has to be followed by `-n' and the weight file names for dim>4");
}

/*  ==========             ALL  IP  WEIGHTS  in  d <= 4           ==========  */
#define  WDIM      300000
#define  INFO      0

#define  lcm(a,b)  ((a)*(b)/NNgcd((a),(b)))
typedef struct {int n[W_Nmax+1];} weights;
typedef Rat ratmat[W_Nmax][W_Nmax];
typedef Rat ratvec[W_Nmax];
typedef struct { int wnum, N, points[W_Nmax][W_Nmax], nsubsets[W_Nmax-1],
      subsets[W_Nmax-1][10][W_Nmax];  weights wli[WDIM]; }  WSaux;

int  weicomp(weights w1,weights w2,int *_N) 
                                /* w2-w1, i.e. pos for w1<w2,neg for w1>w2  */
{    int i=*_N;
     while((i)&&(w1.n[i]==w2.n[i])) i--;  return w2.n[i]-w1.n[i];
}
void insertat(WSaux *X,weights ww, int position)
{    int i, j;
     for(i=X->wnum-1;i>=position;i--)
       for(j=0;j<=X->N;j++) X->wli[i+1].n[j]=X->wli[i].n[j];
     for(j=0;j<=X->N;j++) X->wli[position].n[j]=ww.n[j];     X->wnum++;
}

void addweight(WSaux *X,weights wn)
{    int i, p, n0, n1, k;
     if (X->wnum>=WDIM) {
       if(X->wnum>WDIM) return; 
       else {X->wnum++; printf("WDIM too small!\n");fflush(0);return;}}
     for(i=0;i<X->N-1;i++) for(p=i+1;p<X->N;p++) if (wn.n[i]>wn.n[p]) 
       {k=wn.n[i]; wn.n[i]=wn.n[p]; wn.n[p]=k;}  /* make n0<=n1<=...<=n# */
     if (X->wnum) {if ((i = weicomp(wn,X->wli[n0=0],&X->N)))
               {if (i>0) {insertat(X,wn,0); return;}} else return;
            if ((i = weicomp(wn,X->wli[n1=X->wnum-1],&X->N)))
               {if (i<0) {insertat(X,wn,X->wnum); return;}} else return;
      while(n1>n0+1) {p=(n0+n1)/2; i=weicomp(wn,X->wli[p],&X->N);
                  if(i) {if(i>0) n1=p; else n0=p;} else return;}
      insertat(X,wn,n1);}
     else insertat(X,wn,0);
}
int checkwrite(WSaux *X,weights ws){
  int i;
  for (i=0;i<X->N;i++) if (!ws.n[i]) return 0;
  addweight(X,ws);
  return 1;
}
weights testweisys(WSaux *X,int npoints){
  weights tws;
  ratmat rm;
  ratvec newboundwei, boundwei[10], rattws, rs;
  int i, ii, j, k, nboundwei=0, New, rankrm, one[W_Nmax];
  Long minnbw, maxnbw;
  for (k=0;k<X->nsubsets[npoints-2];k++   /* alle 0-systeme */ ){
    rankrm=0;
    for (i=0;i<npoints;i++){
      for (j=0;j<npoints;j++){
        rm[i][j]=rI(X->points[i][X->subsets[npoints-2][k][j]]);}
      rs[i]=rI(1);   }
    for (i=0;i<npoints;i++){
      one[i]=-1;
      for (j=0;(j<npoints)&&(one[i]<0);j++) if (rm[i][j].N) {
        one[i]=j;
        rankrm++;   }
      if (one[i]>=0){
        for (j=one[i]+1;j<npoints;j++) /* normalize i'th line */
          rm[i][j]=rQ(rm[i][j],rm[i][one[i]]);
        rs[i]=rQ(rs[i],rm[i][one[i]]);
        rm[i][one[i]]=rI(1);
        for (ii=0;ii<npoints;ii++) if (ii-i){ /* subtract multiple of i'th
                    line from ii'th line */
          for (j=0;j<npoints;j++) 
            if (j!=one[i]) rm[ii][j]=rD(rm[ii][j],rP(rm[ii][one[i]],rm[i][j]));
          rs[ii]=rD(rs[ii],rP(rm[ii][one[i]],rs[i]));
          rm[ii][one[i]]=rI(0);  }   }    }
    for (j=0;j<X->N;j++) newboundwei[j]=rI(0);
    if (rankrm==npoints) for (i=0;i<npoints;i++) 
      newboundwei[X->subsets[npoints-2][k][one[i]]]=rs[i];
    minnbw=0; maxnbw=0; 
    for (j=0;j<X->N;j++) {
      minnbw=min(newboundwei[j].N,minnbw);
      maxnbw=max(newboundwei[j].N,maxnbw); }
    New=((minnbw>=0)&&(maxnbw>0));
    for (i=0;New&&(i<nboundwei);i++) {
      New=0; 
      for (j=0;j<X->N;j++) if (rD(newboundwei[j],boundwei[i][j]).N) New=1; }
    if (New) {
      for (j=0;j<X->N;j++) boundwei[nboundwei][j]=newboundwei[j]; 
      nboundwei++;    }   }
  for (j=0;j<X->N;j++) {
    rattws[j]=rI(0);
    for (i=0;i<nboundwei;i++) rattws[j]=rS(rattws[j],boundwei[i][j]);
    rattws[j]=rQ(rattws[j],rI(max(nboundwei,1)));}
  tws.n[X->N]=1;
  for (j=0;j<X->N;j++) tws.n[X->N]=lcm(tws.n[X->N],rattws[j].D);
  for (j=0;j<X->N;j++) tws.n[j]=rP(rI(tws.n[X->N]),rattws[j]).N;
  return tws;
}
void createweights(WSaux *X,int npoints){
  int x0, x1, x2, x3, x4, sum, maxx;  weights tws; tws=testweisys(X,npoints);
  if(checkwrite(X,tws)) if (npoints<X->N){
    for (x0=0;x0*tws.n[0]<tws.n[X->N];x0++)
    for (x1=0;x0*tws.n[0]+x1*tws.n[1]<tws.n[X->N];x1++)
    for (x2=0;x0*tws.n[0]+x1*tws.n[1]+x2*tws.n[2]<tws.n[X->N];x2++)
/* #if (N>3) */    for (x3=0; (x3==0) || ((X->N>3) &&
      (x0*tws.n[0]+x1*tws.n[1]+x2*tws.n[2]+x3*tws.n[3]<tws.n[X->N])) ;x3++)
/* #if (N>4) */    for (x4=0; (x4==0) || ((X->N>4) && (x0*tws.n[0]+
      x1*tws.n[1]+x2*tws.n[2]+x3*tws.n[3]+x4*tws.n[4]<tws.n[X->N]));x4++)
    { sum=0; maxx=0;
      X->points[npoints][0]=x0; sum+=x0; maxx=max(maxx,x0);
      X->points[npoints][1]=x1; sum+=x1; maxx=max(maxx,x1);
      X->points[npoints][2]=x2; sum+=x2; maxx=max(maxx,x2);
  if (X->N>3)    
      {X->points[npoints][3]=x3; sum+=x3; maxx=max(maxx,x3);}
  if (X->N>4)    
      {X->points[npoints][4]=x4; sum+=x4; maxx=max(maxx,x4);}
      if ((sum>2)&&(maxx>1)) createweights(X,npoints+1);
      /* if (npoints<3) {printf("%d",npoints); fflush(0);}*/}  }
}
void makesubsets(WSaux *X){
  int i, p0,p1,p2,p3,p4;
  for (i=0;i<X->N-1;i++) X->nsubsets[i]=0;
  for (p0=0;p0<X->N-1;p0++) for (p1=p0+1;p1<X->N;p1++) {
    X->subsets[0][X->nsubsets[0]][0]=p0; 
    X->subsets[0][X->nsubsets[0]][1]=p1;
    X->nsubsets[0]++;
    for (p2=p1+1;p2<X->N;p2++) {
      X->subsets[1][X->nsubsets[1]][0]=p0; 
      X->subsets[1][X->nsubsets[1]][1]=p1;
      X->subsets[1][X->nsubsets[1]][2]=p2;
      X->nsubsets[1]++;
      for (p3=p2+1;p3<X->N;p3++) {
        X->subsets[2][X->nsubsets[2]][0]=p0; 
        X->subsets[2][X->nsubsets[2]][1]=p1;
        X->subsets[2][X->nsubsets[2]][2]=p2;
        X->subsets[2][X->nsubsets[2]][3]=p3;
        X->nsubsets[2]++;
        for (p4=p3+1;p4<X->N;p4++) {
          X->subsets[3][X->nsubsets[3]][0]=p0; 
          X->subsets[3][X->nsubsets[3]][1]=p1;
          X->subsets[3][X->nsubsets[3]][2]=p2;
          X->subsets[3][X->nsubsets[3]][3]=p3;
          X->subsets[3][X->nsubsets[3]][4]=p4;
          X->nsubsets[3]++;  }  }  }  } 
}
void WRITE_Weight(Weight *_W);
void Make_34_Weights(int d, int tFlag)
{    int i, Info=0; WSaux *X = (WSaux *) malloc(sizeof(WSaux)); 
     PolyPointList *P = (PolyPointList *) malloc(sizeof (PolyPointList));
     assert(P!=NULL); assert(X!=NULL); X->wnum=1; X->N=d+1; assert(d<=4); 
     makesubsets(X);    for (i=0;i<X->N;i++) X->points[0][i]=1;
     X->wli[0].n[X->N]=X->N;  for (i=0;i<X->N;i++) X->wli[0].n[i]=1;
  if (X->N>4){ X->points[1][0]=4; for (i=1;i<X->N;i++) X->points[1][i]=0;
    createweights(X,2); if(Info){printf("Did (4,0,0,0,0)\n");fflush(0);} 
    X->points[1][0]=3; X->points[1][1]=1; for(i=2;i<X->N;i++)X->points[1][i]=0;
    createweights(X,2);  if(Info){printf("Did (3,1,0,0,0)\n"); fflush(0);}
    X->points[1][0]=2; X->points[1][1]=2; for(i=2;i<X->N;i++)X->points[1][i]=0;
    createweights(X,2);  if(Info){printf("Did (2,2,0,0,0)\n"); fflush(0);}
    X->points[1][0]=2; X->points[1][1]=1; X->points[1][2]=1;
      for (i=3;i<X->N;i++) X->points[1][i]=0;
    createweights(X,2);  if(Info){printf("Did (2,1,1,0,0)\n"); fflush(0);}}
  if (X->N>3){ X->points[1][0]=3; for (i=1;i<X->N;i++) X->points[1][i]=0;
    createweights(X,2);  if(Info){printf("Did (3,0,0,0,0)\n"); fflush(0);}
    X->points[1][0]=2; X->points[1][1]=1; for(i=2;i<X->N;i++)X->points[1][i]=0;
    createweights(X,2);  if(Info){printf("Did (2,1,0,0,0)\n"); fflush(0);}}         X->points[1][0]=2; for (i=1;i<X->N;i++) X->points[1][i]=0;
  createweights(X,2);  if(Info){printf("Did (2,0,0,0,0)\n"); fflush(0);}
  Info=0; for (i=0;i<X->wnum;i++){ int j; Weight W; VertexNumList V; EqList E;
    W.d=X->wli[i].n[W.N=X->N]; for(j=0;j<X->N;j++)W.w[j]=X->wli[i].n[j]; W.M=0;
    Make_Poly_Points(&W, P); if(Ref_Check(P,&V,&E))
    { int t=Trans_Check(W); char c[5]="  rt"; if(t||!tFlag) {c[3]=(t)?'t':0; 
      if(Info++) puts(""); WRITE_Weight(&W); printf("%s",c);}}
  } fprintf(outFILE,"  #=%d  #cand=%d\n",Info,X->wnum); free(X);
}
/*  ==========       End of  ALL  IP  WEIGHTS  in  d <= 4         ==========  */
/*  ==========          MAKE WEIGHTS d>4:                   ==========  */

void WRITE_Weight(Weight *_W)
{    int n; fprintf(outFILE,"%d ",(int) _W->d);
     for(n=0; n < _W->N; n++) fprintf(outFILE," %d",(int) _W->w[n]); 
}

int IfIpWWrite(Weight *W, PolyPointList *P, int *rFlag, int *tFlag)
{    VertexNumList V; EqList E; int r=1,i=-1; Make_Poly_Points(W, P); 
     if(IP_Check(P,&V,&E)){
       while(r && (++i < E.ne)) if(E.e[i].c != 1) r = 0;
       if(*tFlag && Trans_Check(*W)){
       WRITE_Weight(W); if(r) fprintf(outFILE," r"); 
       fprintf(outFILE,"\n");fflush(stdout); return 1;}
       if(*rFlag && r){Write_Weight(W); fflush(stdout); return 1;} 
       if(!*tFlag && !*rFlag)
       {WRITE_Weight(W); if(r) fprintf(outFILE," r"); 
#if(TRANS_INFO_FOR_IP_WEIGHTS)
        if(Trans_Check(*W)) fprintf(outFILE,"%st", r ? "" : " "); 
#endif
       fprintf(outFILE,"\n"); fflush(stdout); return 1;}
       return 0;
     } 
     else return 0;
}
void Rec_IpWeights(Weight *W, PolyPointList *P, int g, int sum, int *npp, 
      int *nrp, int n, int *rFlag, int *tFlag)
{    int wmax=W->d/(W->N-n+1); wmax=min(wmax,W->w[n+1]); 
     wmax=min(wmax,sum-n);
     if(n) for(W->w[n]=wmax;(n+1)*W->w[n]>=sum;W->w[n]--)
       Rec_IpWeights(W,P,Fgcd(g,W->w[n]),sum-W->w[n],npp,nrp,n-1,rFlag,tFlag);
     else if(1==Fgcd(g,W->w[0]=sum)) {
       (*npp)++;if(IfIpWWrite(W,P,rFlag,tFlag))(*nrp)++;};
}
void MakeIpWeights(int N, int from_d, int to_d, int *rFlag, int *tFlag)
{    int npp=0, nrp=0; Weight W;  
     PolyPointList *P = (PolyPointList *) malloc (sizeof(PolyPointList)); 
     assert((N<=W_Nmax)&&(N<POLY_Dmax+2)); assert(P!=NULL);W.N=N; W.M=0;
     for(W.d=from_d;W.d<=to_d;W.d++)
        for(W.w[N-1]=W.d/2; W.d <= N*W.w[N-1]; W.w[N-1]--)
      Rec_IpWeights(&W,P,Fgcd(W.d,W.w[W.N-1]), W.d-W.w[W.N-1],&npp,&nrp,
                  N-2,rFlag,tFlag);
     if(*rFlag) fprintf(outFILE,"#primepartitions=%d #refpolys=%d\n",npp,nrp);
    if(*tFlag)fprintf(outFILE,"#primepartitions=%d #transpolys=%d\n",npp,nrp);
     if(!*rFlag && !*tFlag)
       fprintf(outFILE,"#primepartitions=%d #IPpolys=%d\n",npp,nrp);
     exit(0);
}
void Make_IP_Weights(int d, int Dmin, int Dmax, int rFlag, int tFlag)
{
  MakeIpWeights(d+1, Dmin, Dmax, &rFlag, &tFlag);
}
      /* ----------   LG/transversal stuff  ------------ */
#define  ALLOWHALF      (1)               /* i.e. trivial LG potentials */
#define  CHAT           (0)               /* 3 ... for positive c_1 */
#define  TWDIM  16384      /* 16384  8192  4096  dimension of weight-buffer */
#define  mod(a,b)  ((a)%(b))

typedef int T_weight[AMBI_Dmax+2];              /* NM::AMBI_Dmax */
typedef struct {int n,d,wnum,jmax;T_weight wei, wli[TWDIM];}      T_aux;

void T_Chon(int,int,int,int,T_aux*); /* i, {-fermat,0=closed,+open}, nmax, g */
void T_Addweight(T_weight,T_aux *X);
int PPT_Check(T_weight nli,T_aux *X);

void Make_Trans_Weights(int n,int dmin,int dmax /*,int rFlag */)
{    int i,j,inc=1; T_aux X; X.n=n; X.wnum=0;
     outFILE=stdout; assert(n<=AMBI_Dmax); X.wei[0]=n;
     if(CHAT){ assert(CHAT==3); if(!(n%2)) {inc++; dmin+=(dmin%2);}}
      for(X.d=dmin;X.d<=dmax;X.d+=inc) 
      { X.wei[n+1]=X.d; X.wnum=0;
      if(ALLOWHALF) X.jmax=X.d/2; else X.jmax=(X.d-1)/2; 
      if(CHAT) T_Chon(1,-1,(X.d*(n- CHAT))/2-n+1,X.d,&X);  
      else T_Chon(1,-1, X.d -n+1,X.d,&X);
        for(i=0;i<X.wnum;i++) 
        {   Weight W; W.N=n; W.d=X.d; for(j=0;j<n;j++) W.w[j]=X.wli[i][j+1];
            W.M=0; if(Trans_Check(W)) Write_Weight(&W);       
        }
        fflush(outFILE);
      }
}
/* T_Chon chooses 0<wei[i]<=nm=d(n-3)/2-n+i-wei[1]-...-wei[i-1] such that a *
 * pointer structur with at most one unresolved pointer(urp) pointing to the*
 * right at a time is not ruled out. If (urp!=0) then the number of the un- *
 * resolved pointer is urp; in addition put all fermats to the left (urp<0).*
 * i has to point at l>=urp                                                 */
/* let j run; check mod(d||(d-n),j); if (upr) check if upr is resolved by j;*/
void T_Chon(int i, int urp, int nm, int g,T_aux *X)
{    int res, j, l=0, ip=i+1, jm=min(nm,X->jmax);
     if (i<X->n) for(j=(i==X->n-1) ? (1+nm-jm) : 1;j<=jm;j++)
            {X->wei[i]=j;/*next step*/
        if(urp<0) {if(X->d%j) {res=1; /* i.e. not ferm; res=0 -> resolved */
                      for(l=1;(l<i)&&res;l++) res=(X->d-X->wei[l])%j;
                      if (res) T_Chon(ip,i,nm-j+1,Fgcd(g,j),X);
                      else T_Chon(ip,0,nm-j+1,Fgcd(g,j),X);}
                 else{if((i==1)||(j>=X->wei[i-1]))
                  T_Chon(ip,-1,nm-j+1,Fgcd(g,j),X);}
        } else  /* continue;} */
        if(X->d%j) {                   /* now there can be no more fermat */
        l=max(urp,1); for(res=1;(l<i)&&res;l++) res=(X->d-X->wei[l])%j;
        if(urp){if(!res){if((X->d-j)%X->wei[urp]) 
                  T_Chon(ip,urp,nm-j+1,Fgcd(g,j),X);
                            else T_Chon(ip,0,nm-j+1,Fgcd(g,j),X);}}
        else {if(res) T_Chon(ip,i,nm-j+1,Fgcd(g,j),X);
              else T_Chon(ip,0,nm-j+1,Fgcd(g,j),X);} }}
     else {if(X->jmax<nm) return;                              /* last step */
        if(urp>0) if((X->d-nm)%X->wei[urp]) return;    /* pointer resolved? */
        if((res=(X->d%nm))||(urp<0)){ 
            for(l=1;(l<i)&&res;l++) res=(X->d-X->wei[l])%nm;
        if(!res){X->wei[i]=nm;if(Fgcd(g,nm)==1)PPT_Check(X->wei,X);}} }
}
/*  ppcheck checks whether the formal poincare polynomial is a polynomial  */
int PPT_Check(T_weight nli,T_aux *X)
{    int i=0, n, t, tt, j,d=X->d; assert(d==nli[nli[0]]);
     for (i=1;i<=nli[0];i++){
        n=1; tt=nli[i];
        for (j=i+1;j<=nli[0];j++) if (!(nli[j]%tt)) n++;
        for (j=1;j<=nli[0];j++) if (!((d-nli[j])%tt)) n--;
        if (n>0) return 0;
        for (t=2;t*t<=nli[i];t++)
           if (!mod(nli[i],t)) {
              n=1;
              for (j=i+1;j<=nli[0];j++) if (!mod(nli[j],t)) n++;
              for (j=1;j<=nli[0];j++) if (!mod(d-nli[j],t)) n--;
              if (n>0) return 0;
              n=1; tt=nli[i]/t;
              for (j=i+1;j<=nli[0];j++) if (!mod(nli[j],tt)) n++;
              for (j=1;j<=nli[0];j++) if (!mod(d-nli[j],tt)) n--;
              if (n>0) return 0;};};
     T_Addweight(nli,X); return 1;
}
int  T_Weicomp(T_weight w1,T_weight w2)/* w2-w1,i.e.pos if w1<w2,neg if w1>w2*/
{    int i=1; while((i<=(*w1)) && (w1[i]==w2[i])) i++; return w2[i]-w1[i];
}
void T_Insertat(T_weight ww, int position,T_aux *X)
{    int i, j;
     for(i=X->wnum-1;i>=position;i--)
         for(j=0;j<X->wli[i][0]+2;j++) X->wli[i+1][j]=X->wli[i][j];
     for(j=0;j<ww[0]+2;j++) X->wli[position][j]=ww[j]; X->wnum++;
}
void T_Addweight(T_weight win,T_aux *X)
{    int i, p, n0, n1; T_weight wn; for(i=0;i<*win+2;i++) wn[i]=win[i];
     for(i=1;i<wn[0];i++) for(p=i+1;p<=wn[0];p++)
             if (wn[i]>wn[p]) swap(&wn[i],&wn[p]);  /* make n0<=n1<=...<=n# */
     if (X->wnum) {if ((i = T_Weicomp(wn,X->wli[n0=0])))
                   {if (i>0) {T_Insertat(wn,0,X); return;}} else return;
                if ((i = T_Weicomp(wn,X->wli[n1=X->wnum-1])))
                   {if (i<0) {T_Insertat(wn,X->wnum,X); return;}} else return;
        while(n1>n0+1) {p=(n0+n1)/2; i=T_Weicomp(wn,X->wli[p]);
        if(i) {if(i>0) n1=p; else n0=p;} else return;}
        T_Insertat(wn,n1,X);}
     else T_Insertat(wn,0,X);
}
/*  ==========      End of MAKE WEIGHTS d>4:          ==========  */




/*  ==========              ALL  CWS  in  d <= 4            ==========  */
typedef struct {int d, w[2];}   wei2; 
typedef struct {int d, w[3];}   wei3; 
typedef struct {int d, w[4];}   wei4;
const wei2 W2={2,{1,1}};
const wei3 W3[ 3]={{3,{1,1,1}},   {4,{1,1,2}},      {6,{1,2,3}}};
const wei4 W4[95]={
{4, {1,1,1,1}},    {5,{1,1,1,2}},   {6,{1,1,2,2}},    {6,{1,1,1,3}},    {7,{1,1,2,3}},
{8, {1,2,2,3}},    {8,{1,1,2,4}},   {9,{1,2,3,3}},    {9,{1,1,3,4}},    {10,{1,2,3,4}},
{10,{1,2,2,5}},    {10,{1,1,3,5}},{11,{1,2,3,5}},     {12,{2,3,3,4}},   {12,{1,3,4,4}},
{12,{2,2,3,5}},    {12,{1,2,4,5}},{12,{1,2,3,6}},     {12,{1,1,4,6}},   {13,{1,3,4,5}},
{14,{2,3,4,5}},    {14,{2,2,3,7}},{14,{1,2,4,7}},     {15,{3,3,4,5}},   {15,{2,3,5,5}},
{15,{1,3,5,6}},    {15,{1,3,4,7}},{15,{1,2,5,7}},     {16,{1,4,5,6}},   {16,{2,3,4,7}},
{16,{1,3,4,8}},    {16,{1,2,5,8}},{17,{2,3,5,7}},     {18,{3,4,5,6}},   {18,{1,4,6,7}},
{18,{2,3,5,8}},    {18,{2,3,4,9}},{18,{1,3,5,9}},     {18,{1,2,6,9}},   {19,{3,4,5,7}},
{20,{2,5,6,7}},    {20,{3,4,5,8}},{20,{2,4,5,9}},     {20,{2,3,5,10}},{20,{1,4,5,10}},
{21,{3,5,6,7}},    {21,{1,5,7,8}},{21,{2,3,7,9}},     {21,{1,3,7,10}},{22,{2,4,5,11}},
{22,{1,4,6,11}}, {22,{1,3,7,11}},{24,{3,6,7,8}},{24,{4,5,6,9}},   {24,{1,6,8,9}},
{24,{3,4,7,10}},{24,{2,3,8,11}},{24,{3,4,5,12}},{24,{2,3,7,12}},{24,{1,3,8,12}},
{25,{4,5,7,9}},   {26,{2,5,6,13}},{26,{1,5,7,13}},{26,{2,3,8,13}},{27,{5,6,7,9}},
{27,{2,5,9,11}},{28,{4,6,7,11}},{28,{3,4,7,14}},{28,{1,4,9,14}},{30,{5,6,8,11}},
{30,{3,4,10,13}}, {30,{4,5,6,15}}, {30,{2,6,7,15}}, {30,{1,6,8,15}},
{30,{2,3,10,15}}, {30,{1,4,10,15}},{32,{4,5,7,16}}, {32,{2,5,9,16}},
{33,{3,5,11,14}}, {34,{4,6,7,17}}, {34,{3,4,10,17}},{36,{7,8,9,12}},
{36,{3,4,11,18}}, {36,{1,5,12,18}},{38,{5,6,8,19}}, {38,{3,5,11,19}},
{40,{5,7,8,20}},  {42,{3,4,14,21}},{42,{2,5,14,21}},{42,{1,6,14,21}},
{44,{4,5,13,22}}, {48,{3,5,16,24}},{50,{7,8,10,25}},{54,{4,5,18,27}},
{66,{5,6,22,33}}};

void MakeSelections(FILE *, FILE *, int);
void Make2CWS(FILE *, FILE *, int, int);
void RW_TO_CWS(CWS *, Weight *, int, int, int, int);
void W_TO_CWS(CWS *, Weight *, int, int, int, int);
void PRINT_CWS(CWS *);
void Make_111_CWS(FILE **, int *);
void Make_nno_CWS(FILE **, int, int);

void STtmp(FILE *w2FILE, FILE *w3FILE, FILE *w4FILE)
{
  int i,j;

  for(i=0; i<95; i++){fprintf(w4FILE, "%d  ",W4[i].d);
    for(j=0; j<4; j++) fprintf(w4FILE, "%d ",W4[i].w[j]);
    fprintf(w4FILE,"\n");
  }
  for(i=0; i<3; i++){fprintf(w3FILE, "%d  ",W3[i].d);
    for(j=0; j<3; j++) fprintf(w3FILE, "%d ",W3[i].w[j]);
    fprintf(w3FILE,"\n");
  }
  fprintf(w2FILE, "%d  ",W2.d); for(j=0;j<2;j++)fprintf(w2FILE, "%d ",W2.w[j]);
  fprintf(w2FILE,"\n");
  rewind(w2FILE); rewind(w3FILE); rewind(w4FILE);
}

void mkold2(char *outfile, FILE *INFILE1, FILE *INFILE2, int u, int ef)
{
  
  FILE *AUXFILE1, *AUXFILE2;
  
  if((AUXFILE1 = tmpfile()) == NULL) 
      Die("Unable to open tmpfile for read/write");
  if((AUXFILE2 = tmpfile()) == NULL) 
      Die("Unable to open tmpfile for read/write");

  MakeSelections(INFILE1, AUXFILE1, u);
  MakeSelections(INFILE2, AUXFILE2, u);

  if(strcmp(outfile,""))
    if((outFILE = fopen(outfile, "w")) == NULL){
      printf("\nUnable to open file %s for write\n",outfile); exit(0);}
 
  Make2CWS(AUXFILE1, AUXFILE2, u, ef);
  fclose(AUXFILE1); fclose(AUXFILE2);
  if(strcmp(outfile,""))
    fclose(outFILE);
}

void mk2xxx(char *outfile, int n)
{
  int i,j,d=0;
  CWS CW; CW.nz=0;

  CW.N=2*n; CW.nw=n; 
  for(i=0; i<n; i++){
    CW.d[i]=2; 
    for(j=0; j<CW.N; j++){
      if((j==d) || (j==d+1)) CW.W[i][j]=1;
      else CW.W[i][j]=0;
    }
    d=d+2; 
  }
  if(strcmp(outfile,""))
    if((outFILE = fopen(outfile, "w")) == NULL){
      printf("\nUnable to open file %s for write\n",outfile); exit(0);}
  PRINT_CWS(&CW);
  if(strcmp(outfile,""))
    fclose(outFILE);
}

void mk3u3u3(char *outfile, FILE *INFILE)
{
 
  FILE *AUXFILE[3] = {NULL};
  int i, u = 1, eq[2]; 

  for(i = 0; i < 3; i++){
    if((AUXFILE[i] = tmpfile()) == NULL) 
      Die("Unable to open tmpfile for read/write");
    MakeSelections(INFILE, AUXFILE[i], u);
  }
  if(strcmp(outfile,""))
    if((outFILE = fopen(outfile, "w")) == NULL){
      printf("\nUnable to open file %s for write\n",outfile); exit(0);}
 
  eq[0] = eq[1] = 1;
  Make_111_CWS(AUXFILE, eq);
  for(i = 0; i < 3; i++) fclose(AUXFILE[i]); 
  if(strcmp(outfile,""))
    fclose(outFILE);
}

void mkold_nno(char *outfile, FILE *INFILE1, FILE *INFILE2, FILE *INFILE3, 
      int u, int eq)
{
  FILE *AUXFILE[3]; 
  
  if((AUXFILE[0]=tmpfile())==NULL)Die("Unable to open tmpfile for read/write");
  if((AUXFILE[1]=tmpfile())==NULL)Die("Unable to open tmpfile for read/write");

  MakeSelections(INFILE1, AUXFILE[0], u);
  MakeSelections(INFILE2, AUXFILE[1], u);

  AUXFILE[2] = INFILE3;
  if(strcmp(outfile,""))
    if((outFILE = fopen(outfile, "w")) == NULL){
      printf("\nUnable to open file %s for write\n",outfile); exit(0);}
  
  Make_nno_CWS(AUXFILE, u, eq);
  fclose(AUXFILE[0]); fclose(AUXFILE[1]); 
  if(strcmp(outfile,""))
    fclose(outFILE);
}

void Make_34_CWS(int d)
{    
  FILE *w2FILE, *w3FILE, *w4FILE; 
  int u, ef;
  char *outfile;

  assert(d<=4); /*puts("Implement Make_34_CWS");*/
  if((w2FILE = tmpfile())==NULL) Die("Unable to open tmpfile for read/write");
  if((w3FILE = tmpfile())==NULL) Die("Unable to open tmpfile for read/write");
  if((w4FILE = tmpfile())==NULL) Die("Unable to open tmpfile for read/write");
  
  STtmp(w2FILE, w3FILE, w4FILE);
  
  if(d == 3){outfile="";
  /*outfile="3u3"*/; u = 1; ef = 1; mkold2(outfile, w3FILE, w3FILE, u, ef);
  /*outfile="2x3"*/; u = 0; ef = 0; mkold2(outfile, w2FILE, w3FILE, u, ef);
  /*outfile="2x2x2"*/; mk2xxx(outfile, atoi("3"));
  }
  if(d == 4){ outfile="";
    /*outfile="4uu4";*/ u = 2; ef = 1; mkold2(outfile, w4FILE, w4FILE, u, ef);
    /*outfile="3u4"; */ u = 1; ef = 0; mkold2(outfile, w3FILE, w4FILE, u, ef);
    /*outfile="3x3"; */ u = 0; ef = 1; mkold2(outfile, w3FILE, w3FILE, u, ef);
    /*outfile="2x4"; */ u = 0; ef = 0; mkold2(outfile, w2FILE, w4FILE, u, ef);
    /*outfile="3u3x2";*/ u=1;ef=1;mkold_nno(outfile,w3FILE,w3FILE,w2FILE,u,ef);
    /*outfile="3x2x2";*/ u=0;ef=0;mkold_nno(outfile,w3FILE,w2FILE,w2FILE,u,ef);
    /*outfile="3u3u3";*/            mk3u3u3(outfile, w3FILE);
    /*outfile="2x2x2x2";*/          mk2xxx(outfile, atoi("4"));
  }
  fclose(w2FILE);fclose(w3FILE);fclose(w4FILE);
}
/*  ==========          End of  ALL  CWS  in  d <= 4        ==========  */

/*  ==========            MAKE CWS d>4                ==========  */

void Print_CWS_Zinfo(CWS *CW);
void Print_CWS(CWS *_W)
{
  int i, j;
  
  for (i = 0; i < _W->nw; i++) {
    fprintf(outFILE, "%d ", (int) _W->d[i]);
    for (j = 0; j < _W->N; j++)
      fprintf(outFILE, "%d ", (int) _W->W[i][j]);
    if (i + 1 < _W->nw)
      fprintf(outFILE, " ");
  }  Print_CWS_Zinfo(_W);
}

void print_W(Weight *_s, Weight *_W, FILE *auxFILE){

  int i,j=0;

  fprintf(auxFILE, "%d ", (int) _W->d);
  for(i = 0; i < _s->N; i++)
    fprintf(auxFILE, "%d ", (int) _W->w[_s->w[i]]);
  for(i = 0; i < _W->N; i++)
    if(i != _s->w[j])
      fprintf(auxFILE, "%d ", (int) _W->w[i]);
    else
      if(j < (_s->N - 1))
        j++;
  fprintf(auxFILE, "\n");
}

void next_n(Weight *_s, Weight *_W, int *_n, FILE *auxFILE){
  int i;

  if(_s->N == *_n)
    print_W(_s, _W, auxFILE);
  else
    if(_s->w[_s->N - 1] != (_W->N - 1)){
      if (_W->w[_s->w[_s->N - 1] + 1] == _W->w[_s->w[_s->N - 1]]){
      _s->w[_s->N] = _s->w[_s->N - 1] + 1;
      _s->N ++;
      next_n(_s, _W, _n, auxFILE);
      _s->N --;
      }
      for(i = _s->w[_s->N - 1] + 1; i < _W->N; i++){
      if(_W->w[i] > _W->w[i - 1]){
        _s->w[_s->N] =  i;
        _s->N ++;
        next_n(_s, _W, _n, auxFILE);
        _s->N --;
      }
      }
    }
}

void Select_n_of_W(Weight *_W, int n, FILE *auxFILE){

  int i;
  Weight s;
  
  if(n == 0){
    fprintf(auxFILE, "%d ", (int) _W->d);
    for(i = 0; i < _W->N; i++)
      fprintf(auxFILE, "%d ", (int) _W->w[i]);
    fprintf(auxFILE, "\n");
  }
  else{
    for(i = 1; i < _W->N; i++)
      if(_W->w[i] < _W->w[i-1])
      Die("Weights must be sortet: W_1 <= W_2 <= .... <=W_N!");
    s.N = 1; s.d = _W->d;
    s.w[0] = 0;
    next_n(&s, _W, &n, auxFILE);
    for(i = 1; i < _W->N; i++)
      if(_W->w[i] > _W->w[i-1]){
      s.w[0] = i;
      next_n(&s, _W, &n, auxFILE);
      }
  }
}

void PRINT_CWS(CWS *CW){
#if (!Only_IP_CWS)
  {
    Print_CWS(CW);
    fprintf(outFILE,"\n");
  }
#else
  {
    PolyPointList *P, *DP; EqList E; VertexNumList V;
    P = (PolyPointList *) malloc(sizeof(PolyPointList));
    if (P == NULL) Die("Unable to allocate space for P");
    DP = (PolyPointList *) malloc(sizeof(PolyPointList));
    if (DP == NULL) Die("Unable to allocate space for DP");
    Make_CWS_Points(CW, P);
    if (IP_Check(P,&V,&E)){
      int r=1,i=-1;
      Print_CWS(CW);
      while(r && (++i < E.ne))
      if(E.e[i].c != 1) r = 0;
      Make_Dual_Poly(P, &V, &E, DP);
      fprintf(outFILE," M:%d %d",P->np, V.nv);
      if(r) fprintf(outFILE," N:%d %d",DP->np, E.ne);
      else  fprintf(outFILE," F:%d N:%d", E.ne,DP->np);
      assert(IP_Check(DP,&V,&E));
      fprintf(outFILE,"\n");
    }
    free(DP); free(P); 
  }
#endif  
}

void W_TO_CWS(CWS *CW, Weight *_W, int Nf, int Nr, int Nb, int u)
{  /* d (Nf x 0) w_1 ... w_u  (Nb x 0)  w_(u+1) ... w_N  (Nr x 0) */
  int i, j, d = 0;
    CW->d[CW->nw]=_W->d; CW->N=(Nf+Nb+Nr+_W->N);
  if(CW->N > AMBI_Dmax) Die("increase AMBI_Dmax !");
  for(i=0; i < Nf; i++) CW->W[CW->nw][i]=0;
  for(i=0; i < _W->N; i++){
    CW->W[CW->nw][i+Nf+d]=_W->w[i];
    if(u && (i == (u - 1))){ 
      for(j=0; j < Nb; j++)
      CW->W[CW->nw][i+j+Nf+1]=0;
      d=Nb;
    }
  }
  for(i=0; i < Nr; i++) CW->W[CW->nw][Nf+d+_W->N+i]=0;
  CW->nw++; CW->nz=0;
}

void RW_TO_CWS(CWS *CW, Weight *_W, int Nf, int Nr, int Nb, int u)
{  /* d (Nf x 0) w_N ... w_(u+1)  (Nb x 0)  w_u ... w_1  (Nr x 0) */
  int i, j, d = 0;  
  CW->d[CW->nw]=_W->d; CW->N=(Nf+Nb+Nr+_W->N);
  if(CW->N > AMBI_Dmax) Die("increase AMBI_Dmax !");
  for(i=0; i < Nf; i++) CW->W[CW->nw][i]=0;
  for(i=0; i < _W->N; i++){
    CW->W[CW->nw][i+Nf+d]=_W->w[_W->N - i - 1];
    if(u && (i == (u - 1))){ 
      for(j=0; j < Nb; j++)
        CW->W[CW->nw][i+j+Nf+1]=0;
      d=Nb;
    }
  }
  for(i=0; i < Nr; i++) CW->W[CW->nw][Nf+d+_W->N+i]=0;
  CW->nw++; CW->nz=0;
}

void SWAP(Long *_A, Long *_B)
{
  Long C;
  C = *_A; *_A = *_B; *_B = C;
}

void scan_dim(int nF, char *infile[], int D[])
{
  int i, j = 0;
  FILE *INfile[NFmax];
  Weight W;

  for(i = 0; i < nF; i++){
    if((INfile[i] = fopen(infile[i], "r"))== NULL){
      printf("\nUnable to open file %s for read\n",infile[i]);exit(0);}
    j = 0; while(READ_Weight(&W, INfile[i]))if (j++) break;
    D[i] = W.N - 1;
    fclose(INfile[i]);
  }
}

int Wcomp(Weight *_W1, Weight *_W2){

  int i,j;
  int A[POLY_Dmax+1], B[POLY_Dmax+1];

  if (_W1->N > (POLY_Dmax+1)) Die("increase POLY_Dmax!");
  if (_W1->N != _W2->N) Die("N1 != N2 in Wcomp!");
      
  for(i = 0; i < _W1->N; i++){
    A[i] = i;
    B[i] = i;
  }
  for(i = 0; i < _W1->N - 1; ++i)
    for(j = _W1->N - 1; j > i; --j){
      if(_W1->w[A[j-1]] >  _W1->w[A[j]])
        swap(&A[j-1], &A[j]);
      if(_W2->w[B[j-1]] >  _W2->w[B[j]])
        swap(&B[j-1], &B[j]);
    }
  for(i = 0; i < _W1->N; i++){
    if(_W1->w[A[i]] >  _W2->w[B[i]])
      return 1;
    if(_W1->w[A[i]] <  _W2->w[B[i]])
      return -1;
  }
  return 0;
}

void Make_nno_CWS(FILE *AUXFILE[], int u, int ef)
{
  int n=0, l[2]; 
  Weight W[3];
  CWS CW;

  l[0]=0; while (READ_Weight(&W[0], AUXFILE[0])){
    l[1]=0; l[0]++; while (READ_Weight(&W[1], AUXFILE[1])){
      l[1]++; while (READ_Weight(&W[2], AUXFILE[2])){
      if((l[0] <= l[1]) || !ef){ CW.nw=0;
        RW_TO_CWS(&CW, &W[0], n, (W[1].N + W[2].N - u), n, n);
        W_TO_CWS(&CW, &W[1], (W[0].N - u), W[2].N, n, n);
        W_TO_CWS(&CW, &W[2], (W[0].N + W[1].N - u), n, n, n);
        PRINT_CWS(&CW);
        if(u == 2)
          if((W[0].w[0] != W[0].w[1]) && (W[1].w[0] != W[1].w[1])){
            SWAP(&W[1].w[0], &W[1].w[1]); CW.nw=0;
            RW_TO_CWS(&CW, &W[0], n, (W[1].N + W[2].N - u), n, n);
            W_TO_CWS(&CW, &W[1], (W[0].N - u), W[2].N, n, n);
            W_TO_CWS(&CW, &W[2], (W[0].N + W[1].N), n, n, n);
            PRINT_CWS(&CW);
          }
      }
      }
      rewind(AUXFILE[2]);
    }
    rewind(AUXFILE[1]);
  }
  rewind(AUXFILE[0]); 
}

void Make_111_CWS(FILE *AUXFILE[], int ef[])
{
  int u = 1, n=0, l[3]; 
  Weight W[3]; 
  CWS CW;

  l[0]=0; while (READ_Weight(&W[0], AUXFILE[0])){
    l[1]=0; l[0]++; while (READ_Weight(&W[1], AUXFILE[1])){
      l[2]=0; l[1]++; while (READ_Weight(&W[2], AUXFILE[2])){
      l[2]++; 
      if(((l[0] <= l[1]) || !ef[0]) && ((l[1] <= l[2]) || !ef[1])){
        CW.nw=0;
        RW_TO_CWS(&CW, &W[0], n, (W[1].N + W[2].N - 2*u), n, n);
        W_TO_CWS(&CW, &W[1], (W[0].N - u), (W[2].N - u), n, n);
        W_TO_CWS(&CW, &W[2], (W[0].N - u), n, (W[1].N - u), u);
        PRINT_CWS(&CW);
      }
      }
      rewind(AUXFILE[2]);
    }
    rewind(AUXFILE[1]);
  }
  rewind(AUXFILE[0]);
}

void Make_222_CWS(FILE *AUXFILE[], int ef[])
{
  int u = 2, n=0, l[3]; 
  Weight W[3]; 
  CWS CW;

  l[0]=0; while (READ_Weight(&W[0], AUXFILE[0])){
    l[1]=0; l[0]++; while (READ_Weight(&W[1], AUXFILE[1])){
      l[2]=0; l[1]++; while (READ_Weight(&W[2], AUXFILE[2])){
      l[2]++; 
      if(((l[0] <= l[1]) || !ef[0]) && ((l[1] <= l[2]) || !ef[1])){
        CW.nw=0;
        RW_TO_CWS(&CW, &W[0], n, (W[1].N + W[2].N - 2*u), n, n);
        W_TO_CWS(&CW, &W[1], (W[0].N - u), (W[2].N - u), n, n);
        W_TO_CWS(&CW, &W[2], (W[0].N - u), n, (W[1].N - u), u);
        PRINT_CWS(&CW);
        if(W[0].w[0] == W[0].w[1]){
          if((W[1].w[0] != W[1].w[1]) && (W[2].w[0] != W[2].w[1])){
            SWAP(&W[1].w[0], &W[1].w[1]); CW.nw=0;
            RW_TO_CWS(&CW, &W[0], n, (W[1].N + W[2].N - 2*u), n, n);
            W_TO_CWS(&CW, &W[1], (W[0].N - u), (W[2].N - u), n, n);
            W_TO_CWS(&CW, &W[2], (W[0].N - u), n, (W[1].N - u), u);
            PRINT_CWS(&CW);
          }
        }
        else{
          if((W[1].w[0] != W[1].w[1]) || (W[2].w[0] != W[2].w[1])){
            SWAP(&W[0].w[0], &W[0].w[1]); CW.nw=0;
            RW_TO_CWS(&CW, &W[0], n, (W[1].N + W[2].N - 2*u), n, n);
            W_TO_CWS(&CW, &W[1], (W[0].N - u), (W[2].N - u), n, n);
            W_TO_CWS(&CW, &W[2], (W[0].N - u), n, (W[1].N - u), u);
            PRINT_CWS(&CW);
          }
          if((W[1].w[0] != W[1].w[1]) && (W[2].w[0] != W[2].w[1])){
            SWAP(&W[1].w[0], &W[1].w[1]); CW.nw=0;
            RW_TO_CWS(&CW, &W[0], n, (W[1].N + W[2].N - 2*u), n, n);
            W_TO_CWS(&CW, &W[1], (W[0].N - u), (W[2].N - u), n, n);
            W_TO_CWS(&CW, &W[2], (W[0].N - u), n, (W[1].N - u), u);
            PRINT_CWS(&CW);
          }
        }
      }
      }
      rewind(AUXFILE[2]);
    }
    rewind(AUXFILE[1]);
  }
  rewind(AUXFILE[0]);
}

void Make_221_CWS(FILE *AUXFILE[], int ef)
{
  int n=0, U = 2, u = 1, i=0, l[2]; 
  Weight W[3]; 
  CWS CW;

  l[0]=0; while (READ_Weight(&W[0], AUXFILE[0])){
    l[1]=0; l[0]++; while (READ_Weight(&W[1], AUXFILE[1])){
      l[1]++; while (READ_Weight(&W[2], AUXFILE[2])){
      if((l[0] <= l[1]) || !ef){
        i=0; CW.nw=0;
        RW_TO_CWS(&CW, &W[0], n, (W[1].N + W[2].N - U -u ), n, n);
        W_TO_CWS(&CW, &W[1], (W[0].N - U), (W[2].N - u ), n, n);
        W_TO_CWS(&CW, &W[2], (W[0].N - U + i), n, (W[1].N - U + 1 - i), u);
        PRINT_CWS(&CW);
        if((W[0].w[0] != W[0].w[1]) && (W[1].w[0] != W[1].w[1])){
          SWAP(&W[1].w[0], &W[1].w[1]); CW.nw=0;
          RW_TO_CWS(&CW, &W[0], n, (W[1].N + W[2].N - U -u ), n, n);
          W_TO_CWS(&CW, &W[1], (W[0].N - U), (W[2].N - u ), n, n);
          W_TO_CWS(&CW, &W[2], (W[0].N - U + i), n, (W[1].N - U + 1 - i), u);
          PRINT_CWS(&CW);
        }
        if((W[0].w[0] != W[0].w[1]) || (W[1].w[0] != W[1].w[1])){
          i=1; CW.nw=0;
          RW_TO_CWS(&CW, &W[0], n, (W[1].N + W[2].N - U -u ), n, n);
          W_TO_CWS(&CW, &W[1], (W[0].N - U), (W[2].N - u ), n, n);
          W_TO_CWS(&CW, &W[2], (W[0].N - U + i), n, (W[1].N - U + 1 - i), u);
          PRINT_CWS(&CW);
          if((W[0].w[0] != W[0].w[1]) && (W[1].w[0] != W[1].w[1])){
            SWAP(&W[1].w[0], &W[1].w[1]); CW.nw=0;
            RW_TO_CWS(&CW, &W[0], n, (W[1].N + W[2].N - U -u ), n, n);
            W_TO_CWS(&CW, &W[1], (W[0].N - U), (W[2].N - u ), n, n);
            W_TO_CWS(&CW, &W[2], (W[0].N - U + i), n,(W[1].N - U + 1 - i),u);
            PRINT_CWS(&CW);
          }
        }
      }
      }
      rewind(AUXFILE[2]);
    }
    rewind(AUXFILE[1]);
  }
  rewind(AUXFILE[0]); 
}

void Make_211_CWS(FILE *AUXFILE[], int ef)
{
  int u = 1, n=0, l[3]; 
  Weight W[3]; 
  CWS CW;

  while (READ_Weight(&W[0], AUXFILE[0])){
    l[1]=0; while (READ_Weight(&W[1], AUXFILE[1])){
      l[2]=0; l[1]++; while (READ_Weight(&W[2], AUXFILE[2])){
      l[2]++; 
      if(((l[1] <= l[2]) || !ef)){ CW.nw=0;
        RW_TO_CWS(&CW, &W[0], n, (W[1].N + W[2].N - 2*u), n, n);
        W_TO_CWS(&CW, &W[1], (W[0].N - u), (W[2].N - u), n, n);
        W_TO_CWS(&CW, &W[2], (W[0].N - 2*u), n, W[1].N, u);
        PRINT_CWS(&CW);
        if((W[0].w[0] != W[0].w[1]) && (W[1].w[0] !=  W[2].w[0])){
          SWAP(&W[0].w[0], &W[0].w[1]); CW.nw=0;
          RW_TO_CWS(&CW, &W[0], n, (W[1].N + W[2].N - 2*u), n, n);
          W_TO_CWS(&CW, &W[1], (W[0].N - u), (W[2].N - u), n, n);
          W_TO_CWS(&CW, &W[2], (W[0].N - 2*u), n, W[1].N, u);
          PRINT_CWS(&CW);
        }
      }
      }
      rewind(AUXFILE[2]);
    }
    rewind(AUXFILE[1]);
  }
  rewind(AUXFILE[0]);
}

void Make2CWS(FILE *AUXFILE1, FILE *AUXFILE2, int u, int ef){

  Weight W1, W2;
  int n=0, l[2]; 
  CWS CW;

  if (u > 2) Die("for u > 2 no support !");
    l[0]=0;while (READ_Weight(&W1, AUXFILE1)){
      l[1]=0; l[0]++;
      while (READ_Weight(&W2, AUXFILE2)){
      l[1]++;
      if((l[1] >= l[0]) || !ef){CW.nw=0;
        RW_TO_CWS(&CW, &W1, n, (W2.N-u), n, n);
        W_TO_CWS(&CW, &W2, (W1.N - u), n, n, n);
        PRINT_CWS(&CW);
        if(u == 2)
          if((W1.w[0] != W1.w[1]) && (W2.w[0] != W2.w[1])){
            SWAP(&W2.w[0], &W2.w[1]); CW.nw=0;
            RW_TO_CWS(&CW, &W1, n, (W2.N - u), n, n);
            W_TO_CWS(&CW, &W2, (W1.N - u), n, n, n);
            PRINT_CWS(&CW);
          }
      }
      }
      rewind(AUXFILE2);
    }
    rewind(AUXFILE1);
}

void MakeSelections(FILE *INFILE, FILE *AUXFILE, int u)
{
  Weight W;

  while (READ_Weight(&W, INFILE))
    Select_n_of_W(&W, u, AUXFILE);
  rewind(AUXFILE); rewind(INFILE);
}

void PrintCWSTypes()
{
  const char B[]="         ";
  printf("\nThe following types are available:\n\n");
  printf("#infiles = 2 (need no -t option):\n");
  printf("%s-c# -n2 [intile1] [infile2] (-t 0 0)\n",B);
  printf("%s-c# -n2 [intile1] [infile2] (-t 1 1)\n",B);
  printf("%s-c# -n2 [intile1] [infile2] (-t 2 2)\n",B);
  printf("#infiles = 3:\n");
  printf("%s-c# -n3 [intile1] [infile2] [infile3] -t n n 0\n",B);
  printf("%s-c# -n3 [intile1] [infile2] [infile3] -t 1 1 1\n",B);
  printf("%s-c# -n3 [intile1] [infile2] [infile3] -t 2 1 1\n",B);
  printf("%s-c# -n3 [intile1] [infile2] [infile3] -t 2 2 1\n",B);
  printf("%s-c# -n3 [intile1] [infile2] [infile3] -t 2 2 2\n",B);
  exit(0);
}

void Make_IP_CWS(int narg, char* fn[])
{
  FILE *INFILE[NFmax] = {NULL}, *AUXFILE[NFmax] = {NULL};
  char *infile[NFmax] = {NULL}, *outfile = NULL, *a;
  int n = 0, d = 0, u = -1, nF = 0, i, D[NFmax];
  CWS_type t;

  while (narg > ++n) {
    if(fn[n][0] != '-') 
      break;
    if(fn[n][1] == 'c'){
      if((fn[n][2]==0) && (narg>n+1)) a=fn[++n]; else a=&fn[n][2];
      if(!IsDigit(*a)) Die("after -c there must be a digit!"); d = atoi(a);
    }
    if(fn[n][1] == 'n'){ 
      if((fn[n][2]==0) && (narg>n+1)) a=fn[++n]; else a=&fn[n][2];
      if(!IsDigit(*a)) Die("after -n there must be a digit!"); nF = atoi(a);
      n++; break;
    }
    if(fn[n][1] == 't') PrintCWSTypes();
  }
  for(i = 0; i < nF; i++){
    if((n >= narg)||(fn[n][0] == '-'))
      {printf("#infiles = %d < %d!\n",i,nF); exit(0);}
    infile[i] = fn[n];
    n++;
  }
  if((narg > n) && (fn[n][0] != '-')) Die("too many infiles!"); n--; t.nu=0;
  while (narg > ++n) {
    if(fn[n][0] != '-') break;
    if(fn[n][1] == 't'){ 
      if((fn[n][2]==0) && (narg>n+1)) a=fn[++n]; else a=&fn[n][2];
      if(*a == 0) PrintCWSTypes();  
      if(!IsDigit(*a)) Die("after -t there must be digit(s)!");  
      t.u[0] = atoi(a); t.nu = 1;
      while ((narg > ++n)  && (t.nu < nF)) {
      if(fn[n][0] == '-') break;
      a=fn[n]; 
      if(!IsDigit(*a)) Die("after -t there must be digit(s)!");
      assert(t.nu < NFmax); t.u[t.nu] = atoi(a); t.nu++;
      }
      n--;
      if(narg > ++n) outfile = fn[n]; 
    }
  }
  if(nF == 0) Die("there is no -n#infiles!"); 
  if(d == 0) Die("No dimensoin specified!");
  if(t.nu && (t.nu != nF)) 
    Die("if input is -nN -t k_1,...,k_n then N must be equal to n!");
  if(outfile == NULL) outFILE = stdout;
  else 
    if((outFILE = fopen(outfile, "w")) == NULL){
      printf("\nUnable to open file %s for write\n",fn[n]); exit(0);}
  scan_dim(nF, infile, D);
  for(i = 0; i < nF; i++){
    if((AUXFILE[i] = tmpfile()) == NULL) 
      Die("Unable to open tmpfile to read/write");
    if((INFILE[i] = fopen(infile[i], "r"))==NULL)
      Die("Unable to open infile to read");
  }
  if(nF == 2){
    if(!t.nu) assert((u = D[0] + D[1] -d) >= 0); 
    else{
      if(t.u[0] != t.u[1]) 
      Die("if input is -n2 -t k_1 k_2 then k_1 must be equal to k_2!");
      if(t.u[0] != (D[0] + D[1] - d))
      Die("wrong DIM -cDIM or wrong TYPES -t TYPE1 TYPE2");
      u = t.u[0];
    }
    for(i = 0; i < nF; i++)
      MakeSelections(INFILE[i], AUXFILE[i], u);
    Make2CWS(AUXFILE[0], AUXFILE[1], u, !strcmp(infile[0], infile[1]));
  }
  if(nF == 3){
    if(!t.nu) Die("with nNUMBER and NUMBER>2 I need -t TYPE1 TYPE2 TYPE3!");
    if(((t.u[0] == 1) && (t.u[1] == 1) && (t.u[2] == 1)) ||
       ((t.u[0] == 2) && (t.u[1] == 2) && (t.u[2] == 2))){ int eq[2];
       if((D[0] + D[1] + D[2] -2*t.u[0] - d) != 0)
      Die("wrong DIM -cDIM or wrong TYPES -t TYPE1 TYPE2 TYPE3");
      eq[0] = 0; eq[1] = 0;
      if(!strcmp(infile[0], infile[1])) eq[0] = 1;
      if(!strcmp(infile[1], infile[2])) eq[1] = 1;
      if(((!eq[0]) && (!eq[1])) && (!strcmp(infile[0], infile[2]))){
      eq[0] = 1;
      MakeSelections(INFILE[0], AUXFILE[0], t.u[0]);
      MakeSelections(INFILE[2], AUXFILE[1], t.u[1]);
      MakeSelections(INFILE[1], AUXFILE[2], t.u[2]);
      }
      else
      for(i = 0; i < nF; i++)
        MakeSelections(INFILE[i], AUXFILE[i], t.u[i]);
      if(t.u[0] == 1)
      Make_111_CWS(AUXFILE, eq);
      else 
      Make_222_CWS(AUXFILE, eq);
    }
    else if(((t.u[0] == 2) && (t.u[1] == 1) && (t.u[2] == 1)) ||
          ((t.u[0] == 2) && (t.u[1] == 2) && (t.u[2] == 1))){
      if((D[0] + D[1] + D[2] -(t.u[1] + t.u[2]) - d) != 0)
      Die("wrong DIM -cDIM or wrong TYPES -t TYPE1 TYPE2 TYPE3");
      for(i = 0; i < nF; i++)
      MakeSelections(INFILE[i], AUXFILE[i], t.u[i]);
      if(t.u[1] == 2)
      Make_221_CWS(AUXFILE, !strcmp(infile[1], infile[2]));
      else
      Make_211_CWS(AUXFILE, !strcmp(infile[0], infile[1]));
    }
    else if((t.u[0] == t.u[1]) && (t.u[2] == 0)){
      if((D[0] + D[1] + D[2] - t.u[0] - d) != 0)
      Die("wrong DIM -cDIM or wrong TYPES -t TYPE1 TYPE2 TYPE3");
      for(i = 0; i < nF; i++)
      MakeSelections(INFILE[i], AUXFILE[i], t.u[i]);
      Make_nno_CWS(AUXFILE, t.u[0], !strcmp(infile[0], infile[1]));
    }
    else PrintCWSTypes();
  }
  for(i = 0; i < nF; i++){fclose(INFILE[i]);fclose(AUXFILE[i]);}
}

/*  ==========          POLY DATA:                    ==========  */

void FileRW(char *file, char *m, FILE *rwFILE){

  if((rwFILE = fopen(file, m)) == NULL){
    printf("\n\nUnable to open file %s for %s!\n\n",file,m);
    exit(0);
  }
}

void IP_Poly_Data(int narg, char* fn[])
{
  int r = 1, i, n = 0, p=0, d=0;
  CWS CW;
  PolyPointList *_P, *_DP;
  VertexNumList *_V;
  EqList *_E;

  inFILE=stdin; outFILE=stdout; /*puts("IP_Poly_Data: to be done");*/

  _P = (PolyPointList *) malloc(sizeof(PolyPointList));
  if (_P == NULL) Die("Unable to allocate space for _P");
  _DP = (PolyPointList *) malloc(sizeof(PolyPointList));
  if (_DP == NULL) Die("Unable to allocate space for _DP");
  _V = (VertexNumList *) malloc(sizeof(VertexNumList));
  if (_V == NULL) Die("Unable to alloc space for VertexNumList _V");
  _E = (EqList *) malloc(sizeof(EqList));
  if (_E == NULL) Die("Unable to alloc space for EqList _E");
  
  


  while ((narg > ++n) && (fn[n][0] == '-')){
    if (fn[n][1] == 'i'){
      if (fn[n][2] == 'p') p=1;
      if (fn[n][2] == 'd') d=1;
    }
    if ((fn[n][1] == 'f') || (fn[n][1] == 0))
      inFILE=NULL;
  }
  n--;
  if (narg > ++n){
    if((inFILE = fopen(fn[n], "r")) == NULL){
      printf("\nUnable to open file %s for read\n",fn[n]);
      exit(0);
    }
  }
  if (narg > ++n){
    if((outFILE = fopen(fn[n], "w")) == NULL){
      printf("\nUnable to open file %s for write\n",fn[n]);
      exit(0);
    }
  }
  while (Read_CWS_PP(&CW, _P))
    if (IP_Check(_P,_V,_E)){
      r=1; i=-1;
      while(r && (++i < _E->ne))
      if(_E->e[i].c != 1)
        r = 0;
      Make_Dual_Poly(_P, _V, _E, _DP);
      if((!p) && (!d)){
      Print_CWS(&CW);
      fprintf(outFILE," M:%d %d",_P->np, _V->nv);
      if(r) fprintf(outFILE," N:%d %d",_DP->np, _E->ne);
      else  fprintf(outFILE," F:%d N:%d", _E->ne,_DP->np);
      }
      if(p) Print_PPL(_P,"");
      if(d) Print_PPL(_DP,"");      
      assert(IP_Check(_DP,_V,_E));
      fprintf(outFILE,"\n");
    }
}
/*  ==========          END POLY DATA                 ==========  */

/*  ==========             Convex Hull:                     ==========  */

int Remove_Identical_Points(PolyPointList *);

int ConvHull(PolyPointList *P1, PolyPointList *P2, PolyPointList *P, 
           VertexNumList *V1, int x)
{
  int i, j;
  VertexNumList *V2;
  EqList *E1, *E2;
  
  E1 = (EqList *) malloc(sizeof(EqList));
  if (E1 == NULL) Die("Unable to alloc space for EqList E1");
  E2 = (EqList *) malloc(sizeof(EqList));
  if (E2 == NULL) Die("Unable to alloc space for EqList E2");
  V2 = (VertexNumList *) malloc(sizeof(VertexNumList));
  if (V2 == NULL) Die("Unable to alloc space for VertexNumList V2");

  Find_Equations(P1, V1, E1); Find_Equations(P2, V2, E2);
  
  if((V1->nv+V2->nv) > VERT_Nmax) Die("increase VERT_Nmax!");
  if((V1->nv+V2->nv) > POINT_Nmax) Die("increase POINT_Nmax!");
  if((P1->n+x) > POLY_Dmax) Die("increase POLY_Dmax!");
  if(x<0) Die("if input is -p#, # must be less than dim of first poly");

  P->np=V1->nv;
  for(i=0; i<V1->nv; i++){
    for(j=P1->n;j<(P2->n+x);j++)P->x[i][j]=0;
    for(j=0;j<P1->n;j++) P->x[i][j]=P1->x[V1->v[i]][j];
  }
  for(i=0; i<V2->nv; i++){
    for(j=0;j<x;j++)P->x[P->np][j]=0;
    for(j=0;j<P2->n;j++) P->x[P->np][x+j]=P2->x[V2->v[i]][j];
    P->np++;
  }
  P->n=P2->n+x;
  if(x==0) Remove_Identical_Points(P);
  i=Find_Equations(P, V1, E1);
  Sort_VL(V1);
  free(E1);free(E2);free(V2);
  return i;
}
void Conv(int narg, char* fn[])
{
  FILE *INFILE[2];
  int n=0, x=0, nF=2, i;
  char *infile[2] = {NULL}, *outfile = NULL, *a;
  PolyPointList *P[2], *PP;
  CWS *CW[2];
  VertexNumList *V;

  V = (VertexNumList *) malloc(sizeof(VertexNumList));
  if (V == NULL) Die("Unable to alloc space for VertexNumList V");
  PP = (PolyPointList *) malloc(sizeof(PolyPointList));
  if (PP == NULL) Die("Unable to allocate space for PolyPointList");

  for(i = 0; i < nF; i++){
    P[i] = (PolyPointList *) malloc(sizeof(PolyPointList));
    if (P[i] == NULL) Die("Unable to allocate space for PolyPointList");
    CW[i] = (CWS *) malloc(sizeof(CWS));
    if (CW[i] == NULL) Die("Unable to allocate space for CWS");
  }

  while ((narg > ++n) && (fn[n][0] == '-')){
    if(fn[n][1] == 'p'){
      if(fn[n][2]!=0){ a=&fn[n][2];
      if(!IsDigit(*a)) Die("after -c there must be a digit!"); x=atoi(a);}
    }
  }
  for(i = 0; i < nF; i++){
    if((n >= narg)||(fn[n][0] == '-'))
      {printf("#infiles = %d < %d!\n",i,nF); exit(0);}
    infile[i] = fn[n];
    n++;
  }
  if(narg > n)
    outfile = fn[n];
  for(i = 0; i < nF; i++)
    if((INFILE[i] = fopen(infile[i], "r"))==NULL)
      Die("Unable to open infile to read");
  if(outfile == NULL) outFILE = stdout;
  else 
    if((outFILE = fopen(outfile, "w")) == NULL){
      printf("\nUnable to open file %s for write\n",fn[n]); exit(0);}
  while(READ_CWS_PP(CW[0], P[0], INFILE[0])){
    while(READ_CWS_PP(CW[1], P[1], INFILE[1]))
      if(ConvHull(P[0], P[1], PP, V, (P[0]->n-x))) Print_VL(PP, V, "Vertices of P");
    rewind(INFILE[1]);
  }
  for(i = 0; i < nF; i++){free(P[i]); free(CW[i]);} free(PP);free(V);  
}
/*  ==========          END of Convex Hull            ==========  */

/*uses latte instead of Aux_Complete_Poly for counting points*/
void td_Print_EL(EqList *_E, int *n, int suppress_c, const char *comment){
  int i,j;
  char command[100];
  sprintf(command,"rm zzL.tmp");
  system(command);
  outFILE=fopen("zzL.tmp","w");
  fprintf(outFILE,"%d %d  %s\n",_E->ne,(*n)+1,comment);
  for(i=0;i<_E->ne;i++) {
    if (!suppress_c) fprintf(outFILE,"%d",(int) _E->e[i].c);
    for(j=0;j<*n;j++) fprintf(outFILE," %3d",(int) _E->e[i].a[j]); 
    fprintf(outFILE,"\n");}
    fclose(outFILE);
}
Long NP_use_lat(EqList *_E, PolyPointList *_P)
{
    int tmp;
    char command[100];
    sprintf(command,"count zzL.tmp | grep '*' | awk '{print $7}' > zzL.tmp1");
    
    td_Print_EL(_E,&_P->n,0,"");
    system(command);outFILE=fopen("zzL.tmp1","r");
    while((fscanf(outFILE,"%d",&tmp))!=EOF);fclose(outFILE); return tmp;
    
}
Long L_Point_Count(Weight *W,PolyPointList *P,VertexNumList *V,EqList *E){
  int i,j,d; 
  Long *G[W_Nmax], GM[W_Nmax][VERT_Nmax];
  d=W->N; for(i=0;i<d;i++) G[i]=GM[i]; W_to_GLZ(W->w,&d, G);
  P->n=d-1; P->np=d; for(i=0;i<d-1;i++)for(j=0;j<d;j++)P->x[j][i]=GM[i+1][j]; 
  Find_Equations(P,V,E); return NP_use_lat(E,P);
}

int Read_Weight(Weight *);
Long Poly_Point_Count(PolyPointList *P,VertexNumList *V,EqList *E);

Long W_Point_Count(Weight *W,PolyPointList *P,VertexNumList *V,EqList *E){
  int i,j,d; 
  Long *G[W_Nmax], GM[W_Nmax][VERT_Nmax];
  d=W->N; for(i=0;i<d;i++) G[i]=GM[i]; W_to_GLZ(W->w,&d, G);
  P->n=d-1; P->np=d; for(i=0;i<d-1;i++)for(j=0;j<d;j++)P->x[j][i]=GM[i+1][j]; 
  return Poly_Point_Count(P,V,E); 
  /*  char c[50]="#points="; sprintf(&c[8],"%d",P->np);  
      if(P->np<20) Print_PPL(P,c); else printf("%s\n",c); */
}
void SimplexPointCount(int narg, char* fn[])
{ Weight W; VertexNumList V; EqList *E = (EqList *) malloc(sizeof(EqList));
  int L; PolyPointList *P = (PolyPointList *) malloc(sizeof(PolyPointList));
  if ((E == NULL)||(P == NULL)) Die("Unable to allocate space for E or P");
  assert(narg>1); if(fn[1][2]=='f') inFILE=NULL; W.M=0; L=(fn[1][1]=='L');
  while(Read_Weight(&W)){int n; Long np= L ? L_Point_Count(&W,P,&V,E) :
    W_Point_Count(&W,P,&V,E); 
    if(np<=SIMPLEX_POINT_Nmax){printf("%d",W.d); 
      for(n=0;n<W.N;n++)printf(" %ld",W.w[n]);printf(" N:%ld\n",np);fflush(0);}}
  free(P); free(E); }

Generated by  Doxygen 1.6.0   Back to index