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

mori.c

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

#undef      FIND_CIRCUITS

FILE *inFILE, *outFILE;

void PrintUsage(){
  puts("\n       ************        USAGE OF mori.x        ************");
  puts("Input: \"PolyPointList\" \"#triangulations\" \"TriangList\" with");
  puts("       TriangList: \"#simpices INCI(simp0) INCI(simp1) ...\""); 
  printf("Output: Mori generators (rank=dim of Mori cone");
  puts(" and Stanley-Reisner ideal\n");
  printf("Example:\necho -e \"2 8\\n0 -1 -1  0 1 -1 -1 0\\n");
  printf("0  2 -1 -1  0  0 1 1\\n1\\n7 10011000 10100100 10110000 10000110");
  printf(" 11000010 10001001 11000001\" | mori.x -f\n");      
  exit(0); 
}

typedef struct { int v, d; Long **x;}                 Matrix;/* Line[v][d] */

void Init_Matrix(Matrix *M,int v, int d){       /* v=#vec, d=dim */
  int i; M->x=(Long **) malloc (v*(sizeof(Long *)+d*sizeof(Long)));
  assert(M->x != NULL); M->d=d; M->v=v;
  M->x[0]=(Long *)(&(M->x[v])); for(i=1;i<v;i++) M->x[i]=&M->x[i-1][d];}

void Free_Matrix(Matrix *M){free(M->x);M->v=M->d=0;}

Long VxV(Long *X,Long *Y,int d){Long z=X[0]*Y[0]; 
  while(--d) z+=X[d]*Y[d]; return z;}

void MxMT(Matrix A,Matrix B,Matrix M){int i,j;    /* M_lin = A_lin x B_col */
  assert((M.v==A.v)&&(M.d==B.v)&&(A.d==B.d));  /* => M.d = B.v    */
  for(i=0;i<M.v;i++)for(j=0;j<M.d;j++)M.x[i][j]=VxV(A.x[i],B.x[j],A.d);}

void Print_LMatrix(Matrix M, char *s){int i,j;
  fprintf(outFILE,"%d %d LV %s\n",M.v,M.d,s); for(i=0;i<M.v;i++){ for(j=0;
    j<M.d;j++) fprintf(outFILE,"%2ld%s",M.x[i][j],(j+1==M.d) ? "\n" : " ");}}
void Print_CMatrix(Matrix M, char *s){int i,j;
  fprintf(outFILE,"%d %d CV %s\n",M.d,M.v,s); for(i=0;i<M.d;i++){ for(j=0;
    j<M.v;j++)fprintf(outFILE,"%2ld%s",M.x[j][i],(j+1==M.v) ? "\n" : " ");}}

Long V_to_GLZ(Long *V,Matrix G){      /* V may contain NULLs, return gcd */
  int i,j,d=G.v,p=0,z=0,*P,*Z = (int*)malloc(d*(2*sizeof(int)+sizeof(Long)));
  Long g,*W; assert(Z!=NULL); P=&Z[d]; W=(Long *) &P[d]; 
  for(i=0;i<d;i++) if(V[i]) {W[p]=V[(P[p]=i)];p++;} else Z[z++]=i; 
  assert(z+p==G.d); 
  if(p>1){ g=W_to_GLZ(W,&p,G.x); if(g<0) for(i=0;i<p;i++) G.x[0][i]*= -1; i=p; 
    while(i--) { int x=P[i]; for(j=p;j<d;j++)G.x[j][x]=0; j=p;
      while(j--) G.x[j][x]=G.x[j][i]; while(0<j--)G.x[j][x]=0; }   
    for(i=0;i<z;i++) for(j=0;j<d;j++) G.x[j][Z[i]] = (d-j==i+1); }
  else { for(j=0;j<d;j++)for(i=0;i<d;i++)G.x[i][j]=(i==j); assert(p); 
    if(P[0]) {G.x[P[0]][P[0]]=G.x[0][0]=0; 
      G.x[0][P[0]]=G.x[P[0]][0]=(V[P[0]]>0) ? 1 : -1;}
    else if(*V<0) G.x[0][0]=-1; g=V[P[0]]; } if(g<0) g=-g;
  for(j=0;j<d;j++) assert(VxV(V,G.x[j],d)==((j==0)*g)); free(Z); 
  return g;}

int  Vdiff_LmR(Long *L,Long *R,int d){ 
  while(d--){Long D=L[d]-R[d]; if(D) return (D>0) ? 1 : -1;} return 0;}

void Aux_G_2_BxG(Matrix G,Matrix B){      /* don't use INV_GLZ::infinite loop */
  int l,c,d=G.d,L=d-B.d; Long *X=B.x[d-1]; assert(L>0);
  for(c=0;c<d;c++) { for(l=L;l<d;l++) { int j; X[l]=0; for(j=L;j<d;j++) 
      X[l]+=B.x[l-L][j-L]*G.x[j][c]; } for(l=L;l<d;l++) G.x[l][c]=X[l];}}
void Aux_MinNonNeg_UT(Matrix M,Matrix G,int c,int r,int d,Long D){
  int i,j; for(i=0;i<r;i++){Long X=VxV(G.x[i],M.x[c],d), R=X/D; 
    if(X-R*D<0)R-=1; for(j=0;j<d;j++)G.x[i][j]-=R*G.x[r][j];}}
int  Make_G_for_GxMT_UT(Matrix M,Matrix G){ /* GxM upper trian, return rank */
  int i,j,r=0,v=M.v,d=M.d; Matrix B; Long *V=(Long*)malloc(d*sizeof(Long));
  assert(G.v==d); assert(G.d==d); assert(V!=NULL); Init_Matrix(&B,d,d); 
  for(i=0;i<d;i++) for(j=0;j<d;j++) G.x[i][j]=(i==j);
  for(i=0;i<v;i++){ int nz=0; 
    for(j=r;j<d;j++) if((V[j]=VxV(G.x[j],M.x[i],d))) nz=1;
    if(nz) {if(r) {Long D; B.d=B.v=d-r; D=V_to_GLZ(&(V[r]),B); 
      Aux_G_2_BxG(G,B); Aux_MinNonNeg_UT(M,G,i,r,d,D);} 
      else V_to_GLZ(V,G); r++;}}
  free(V); Free_Matrix(&B); return r;}          

int  Aux_Is_GLZ(Matrix G, Matrix GI){int i,j,d=G.d;   /* return 0 if not GLZ */
  Matrix M; assert(GI.d==d); Init_Matrix(&M,d,d); 
  for(i=0;i<d;i++) for(j=0;j<d;j++) M.x[j][i]=G.x[i][j];
  Make_G_for_GxMT_UT(M,GI); 
  for(i=0;i<d;i++) for(j=0;j<d;j++) if(VxV(GI.x[i],M.x[j],d)!=(i==j)) 
      {Free_Matrix(&M); return 0;}
  Free_Matrix(&M); return 1;}

void swapL(Long *x,Long *y){Long z=*x;*x=*y;*y=z;}

void Transpose(Matrix M,Matrix MT){int i,j,l=M.v,c=M.d;
  assert((MT.v==c)&&(MT.d==l));
  for(i=0;i<M.v;i++)for(j=0;j<M.d;j++)MT.x[j][i]=M.x[i][j];}

void INV_GLZ(Matrix G,Matrix Ginv){assert(Aux_Is_GLZ(G,Ginv));}

void MoriGen(Matrix T,Long *V){int i; Matrix G; Init_Matrix(&G,T.d,T.d); 
  assert(T.d==T.v+1); assert(T.v==Make_G_for_GxMT_UT(T,G));
  for(i=0;i<T.d;i++) V[i]=G.x[T.v][i]; Free_Matrix(&G);}

void Write_INCI(INCI X,int v);
void FE(char *c){puts("Input format error in ");puts(c);exit(0);}
int  IsDigit(char c) { return (('0'<=c) && (c<='9'));}
int  ReadInt()
{    int n; char c=fgetc(inFILE); if(!IsDigit(c)&&(c!='-')) FE("ReadInt");
     ungetc(c,inFILE); fscanf(inFILE,"%d",&n); 
     while(' '==(c=fgetc(inFILE))); ungetc(c,inFILE); return n;
}
void Read_EOL(){char c;while('\n'!=(c=fgetc(inFILE)))if(feof(inFILE))FE("EOF");}

void INCI_2_VNL(INCI X, VertexNumList *_V, int n){ 
  _V->nv=0; while(!INCI_EQ_0(X)) 
  {n--; if(INCI_M2(X)) _V->v[_V->nv++]=n; X=INCI_D2(X);} Sort_VL(_V);}
void Print_VNL(VertexNumList *V){int i; for(i=0;i<V->nv;i++)
  fprintf(outFILE,"%d ",V->v[i]); fprintf(outFILE,"#=%d\n",V->nv);}

void Print_INCI_list(int nI, INCI *I,int v){int i; printf("INCI[%d]:",nI);
  for(i=0;i<nI;i++) {printf(" ");Write_INCI(I[i],v);} puts("");}

#ifdef      FIND_CIRCUITS                       /* 1<<n == 2^n */
INCI VNL_2_INCI(VertexNumList *V, int n)
{ INCI I=INCI_0(); int i; for(i=0;i<V->nv;i++) I += (1<<(n-1-V->v[i])); 
   /* printf("n=%d I=%ld \n",n,I); Print_VNL(V); INCI_2_VNL(I,V,n);
      Print_VNL(V); exit(0); */                 return I;
  }
#endif

void Print_Mori(PolyPointList *P,int nI, INCI *I){
  int i,j,k,r=0,d=P->n,p=P->np, ngen=0, pli[POLY_Dmax+1],ng0,e0=0,nm=0,nv,np; 
  int m[VERT_Nmax]; Long Z[POLY_Dmax+1]; VertexNumList V; INCI IE[VERT_Nmax];
  PolyPointList *UT = (PolyPointList *) malloc(sizeof(PolyPointList));
  Matrix R,VT,G; EqList *E = (EqList*)malloc(sizeof(EqList)); 
  assert(E!=NULL); assert(UT!=NULL); 
  for(i=1;i<nI;i++)for(j=0;j<i;j++) ngen+=(INCI_abs(INCI_AND(I[i],I[j]))==d); 
  Init_Matrix(&VT,d,d+1);Init_Matrix(&R,ngen,p);Init_Matrix(&G,p,p); ng0=ngen;
  for(i=1;i<nI;i++) for(j=0;j<i;j++) if(INCI_abs(INCI_AND(I[i],I[j]))==d) {
    int a,b,x; INCI_2_VNL(I[j],&V,p); assert(V.nv==d+1); 
    for(k=0;k<=d;k++) pli[k]=V.v[k]; INCI_2_VNL(I[i],&V,p);
    if((V.v[0]==0)&&(pli[0]==0)){k=1; while(pli[k]==V.v[k]) k++; assert(k<=d);}
    else{for(a=0;a<=d;a++)printf("%d ",pli[a]);puts("=Ii");/*adj.simplex pair*/
      for(a=0;a<=d;a++)printf("%d ",V.v[a]);    puts("=Ij ... error");
      Print_INCI_list(nI,I,p); exit(0);}
    if(pli[k]>V.v[k]) pli[0]=V.v[x=k];          /* pli[x-1] < pli[0] < pli[x] */
    else{if(k<d)while(pli[k+1]==V.v[k])if((++k)==d)break; pli[0]=V.v[k];x=k+1;}
//for(a=0;a<=d;a++)printf("%d ",pli[a]);puts("=>pli");
    if(x>1) assert(pli[0]>pli[x-1]); if(x<d) assert(pli[0]<pli[x]);
    for(a=0;a<=d;a++) for(b=0;b<d;b++) VT.x[b][a]=P->x[pli[a]][b]; 
    MoriGen(VT,Z); for(a=0;a<p;a++)R.x[r][a]=0; 
    if(0==Z[0]){ printf("Error: Z[0]==0 for I[%d]&I[%d] !\n",j,i);
      for(a=0;a<=d;a++){printf("r=%d pli[%d]=%d Z[%d]=%ld\n",r,a,pli[a],a,
      Z[a]); fflush(0); assert(pli[a]<p);}exit(0);}
    if(Z[0]>0) for(a=0;a<=d;a++) R.x[r][pli[a]]=Z[a]; 
    else for(a=0;a<=d;a++) R.x[r][pli[a]]=-Z[a];
//for(a=0;a<p;a++)printf("%d ",R.x[r][a]);printf("=Z::R[%d]\n",r);
    for(b=0;b<r;b++) if(!Vdiff_LmR(R.x[r],R.x[b],p)) break; if(b==r) r++; }
  R.v=UT->np=ngen=r; r=Make_G_for_GxMT_UT(R,G); 
  if(r!=P->np-P->n-1) { Matrix GR; R.v=ngen;Print_LMatrix(R,"R");
    Print_LMatrix(G,"GLZ"); Init_Matrix(&GR,R.v,p);
    for(i=0;i<R.v;i++) for(j=0;j<p;j++) GR.x[i][j]=VxV(G.x[j],R.x[i],p);
    Print_LMatrix(GR,"GR"); printf("rank=%d != P.np-P.n-1 !!!\n",r); exit(0);}
//Print_LMatrix(R,"Matrix of all rays"); 
  if(UT->np>=POINT_Nmax){fprintf(outFILE,"need POINT_Nmax>=%d\n",UT->np+1);
    exit(0);} UT->n=r; 
  if(r>POLY_Dmax){fprintf(outFILE,"need POLY_Dmax>=%d\n",UT->n);exit(0);}
  for(i=0;i<ngen;i++) for(j=0;j<r;j++) UT->x[i][j]=VxV(G.x[j],R.x[i],p);
  for(i=0;i<ngen;i++) for(j=r;j<p;j++) assert(0==VxV(G.x[j],R.x[i],p));
  for(i=0;i<UT->n;i++)UT->x[UT->np][i]=0; UT->np++; Find_Equations(UT,&V,E);

  np=UT->np-1; nv=V.nv-1; Sort_VL(&V); 
  if(np!=V.v[nv]) {IDerr(); puts("Suspected INCI data error:");
    Print_INCI_list(nI,I,p); puts("... non-convex triangulation?\n"); exit(0);}
  for(i=0;i<nv;i++) IE[i]=INCI_0();
  for(i=0;i<E->ne;i++) if(Eval_Eq_on_V(&E->e[i],UT->x[np],UT->n)==0){e0++;
    for(j=0;j<nv;j++) IE[j]=INCI_PN(IE[j],
      Eval_Eq_on_V(&E->e[i],UT->x[V.v[j]],r));}/* compute Eq(0)-INCIs for Vs */
  if(e0>VERT_Nmax){fprintf(outFILE,"need VERT_Nmax >= %d\n",e0);exit(0);}
//printf("p=%d nm=%d\n",p,nm);
  for(i=0;i<nv;i++){for(j=0;j<nv;j++) 
    if(INCI_LE(IE[i],IE[j])&&!INCI_EQ(IE[i],IE[j])) break;if(j==nv) m[nm++]=i;}
  assert(nm>=r); fprintf(outFILE,
    "%d mori generators (rank=%d): rays=%d (%d)  #eq=%d (%d)  #v=%d (%d)\n",
    nm,r,ngen,ng0,e0,E->ne,nm,V.nv);
//printf("p=%d nm=%d\n",p,nm);fflush(0);
  for(i=0;i<nm;i++){ int n=V.v[m[i]]; Long s=0; 
    for(j=1;j<p;j++)s+=R.x[n][j]; fprintf(outFILE,"%3ld ",-s);
    for(j=1;j<p;j++)fprintf(outFILE," %2ld",R.x[n][j]); fputs("   I:",outFILE);
    Write_INCI(IE[m[i]],e0); fputs("\n",outFILE);fflush(0);}
  Free_Matrix(&VT); Free_Matrix(&R); Free_Matrix(&G); free(E); free(UT); }

void IFerr(){puts(
  "\n       ************      INPUT FORMAT ERROR        ************");}
void IDerr(){puts(
  "\n       ************       INPUT DATA ERROR         ************");}

INCI Read_INCI(int p)
{    INCI X=INCI_0(); char c; while(' '==(c=fgetc(inFILE))); ungetc(c,inFILE);
     while(IsDigit(c=fgetc(inFILE))) {assert(c<'2'); X=INCI_PN(X,'1'-c); p--;}
     if(p) {IFerr(); printf("Input format error: INCI string too %s\n\n",
      (p>0)? "short" : "long");PrintUsage();} ungetc(c,inFILE); return X;
}
void Write_INCI(INCI X,int v)
{    int i; char c[VERT_Nmax+1]; c[v]=0;
     for(i=0;i<v;i++){c[v-i-1]='0'+INCI_M2(X); X=INCI_D2(X);}
     fprintf(outFILE,"%s",c);
}

void Test_INCI(int *nI, INCI *ILi,int p){INCI X=ILi[0]; VertexNumList V;
  int i; for(i=1;i<*nI;i++)X=INCI_AND(X,ILi[i]);INCI_2_VNL(X,&V,p); 
  if((V.nv!=1)||(V.v[0]!=0)){IDerr();puts("INCI intersection must be 10...0"); 
    exit(0);}X=ILi[0];for(i=1;i<*nI;i++)X=INCI_OR(X,ILi[i]);INCI_2_VNL(X,&V,p);
  if(V.nv!=p) {IDerr();puts("INCI union must be 11...1"); exit(0);}}

void Read_Tri(int p, int *nI, int *nIA, INCI **_I)
{    int i; INCI *I=*_I;assert(0<=(*nIA)); if(*nIA < (*nI = ReadInt())) 
     {      if(0 < *nIA) free(I); *nIA = *nI;               /* realloc I */
      assert(NULL != ( I = (INCI *) malloc(*nIA * sizeof(INCI)) )); *_I=I;
     }      for(i=0;i<*nI;i++) I[i]=Read_INCI(p); Read_EOL(); Test_INCI(nI,I,p);
}
void Print_StanleyReisner(PolyPointList *P,int nI,INCI *I);
void TriList_to_MoriList(PolyPointList *P)
{    int n,nI,nIA=0,p=P->np,Ntri=ReadInt(); static INCI *I=NULL; Read_EOL();
#ifndef     FIND_CIRCUITS
     Print_PPL(P,""); fprintf(outFILE,"%d triangulations:\n",Ntri);fflush(0);
#endif
     for(n=0;n<Ntri;n++) { Read_Tri(p,&nI,&nIA,&I); Print_Mori(P,nI,I);
        Print_StanleyReisner(P,nI,I);} assert(I!=NULL); free(I);
}

/*   =========================================================== */
int main (int narg, char* fn[]){int n=1; CWS *CW=(CWS *) malloc(sizeof(CWS));
  PolyPointList *P = (PolyPointList *) malloc(sizeof(PolyPointList));
  assert((CW!=NULL)&&(P!=NULL)); CW->nw=0;
  { if (narg > 1) if (fn[n][0]=='-') PrintUsage();        /* fn[0]==mori.x */
    if (narg > 1) inFILE=fopen(fn[1],"r"); else inFILE=stdin;
    if (narg > 2) outFILE=fopen(fn[2],"w"); else outFILE=stdout;     
    if(inFILE==NULL) {puts("Cannot open input file"); exit(0);}
    if(outFILE==NULL) {puts("Cannot open output file"); exit(0);}}
// void TestUT(PolyPointList *P);Read_CWS_PP(CW,P);TestUT(P); 
  while(Read_CWS_PP(CW,P)) TriList_to_MoriList(P); free(CW); free(P);return 0;}

void TestUT(PolyPointList *P){
  int i,j; Matrix M,G,T; Init_Matrix(&M,P->np,P->n); 
  for(i=0;i<M.v;i++)for(j=0;j<M.d;j++)M.x[i][j]=P->x[i][j]; 
  Init_Matrix(&G,M.d,M.d); 
  Make_G_for_GxMT_UT(M,G); 
  Print_LMatrix(G,"G"); Print_CMatrix(M,"Pts");
  Init_Matrix(&T,P->n,P->np); MxMT(G,M,T); Print_LMatrix(T,"Triang");
  Free_Matrix(&T);
  {Matrix GI; Init_Matrix(&GI,P->n,P->n); INV_GLZ(G,GI); Free_Matrix(&GI);}
  exit(0);}

#define     SR_GEN_Alloc      (bico*sizeof(INCI))

void Print_StanleyReisner(PolyPointList *P,int nI,INCI *I){long long bico, Abi;
  INCI *IV, *A, *B, *SRG, *G[POLY_Dmax], *M, *N; 
  int i=1,j=P->np/2, m,n, nG[POLY_Dmax+1], g=0,d,v=P->np-1; if(j>P->n)j=P->n+1;
  bico=P->np; while(i<j) {bico*=(P->np-i); bico/=++i;} assert(P->np<VERT_Nmax);
  Abi=bico*sizeof(INCI); A=(INCI *) malloc( Abi ); B=(INCI *) malloc( Abi );
  IV=(INCI *) malloc( v*sizeof(INCI) ); SRG = (INCI *) malloc( SR_GEN_Alloc ),
  assert(SRG!=NULL); assert(IV!=NULL); assert(A!=NULL); assert(B!=NULL);
  IV[0]=INCI_1(); for(i=1;i<v;i++) IV[i]=INCI_PN(IV[i-1],1); M=A; m=0; N=B;
  G[2]=&SRG[g]; nG[2]=0; 
  for(i=1;i<v;i++)for(j=0;j<i;j++) M[m++]=INCI_OR(IV[i],IV[j]);
  for(i=0;i<m;i++){for(j=0;j<nI;j++) if(INCI_LE(M[i],I[j]))break;
    if(j==nI) {SRG[g++]=M[i];nG[2]++;}}                     
  for(d=3;d<=P->n+1;d++){n=0; nG[d]=0; G[d]=&SRG[g];
    for(i=0;i<m;i++){int l=0;
      while(!INCI_LE(IV[l],M[i])) N[n++]=INCI_OR(IV[l++],M[i]);}
    for(i=0;i<n;i++){for(j=0;j<nI;j++) if(INCI_LE(N[i],I[j]))break;
      if(j==nI) {for(j=0;j<g;j++) if(INCI_LE(SRG[j],N[i]))break;
        if(j==g) {SRG[g++]=N[i]; nG[d]++;}}}         assert(n*d==m*(v-d+1));
    INCI *swapI=M; M=N; N=swapI; m=n; }               /* printf("SRI:\n");
             for(i=0;i<m;i++){Write_INCI(M[i],v+1);printf(" ");} puts(""); */
  fprintf(outFILE,"%d Stanley-Reisner generators:\n",g);
  for(i=0;i<g;i++){if(i)fprintf(outFILE," ");Write_INCI(SRG[i],v+1);}
  fputs("\n",outFILE);              free(SRG);free(A);free(B);free(IV);}

Generated by  Doxygen 1.6.0   Back to index