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

Coord.c

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

#undef      TEST_Wbase
#undef  USE_Old_Wbase
#define NO_COORD_IMPROVEMENT        /* switch off weight permutation */

typedef struct {Long x[AMBI_Dmax][AMBI_Dmax]; int n, N;}     CWLatticeBasis;

void Make_CWS_Points(CWS *_C, PolyPointList *_P);
void CWS_to_PermCWS(CWS *Cin, CWS *C, int *pi);

/*  ==========      I/O functions:                    ==========  */

int  IsNextDigit(){    
  char c; c=fgetc(inFILE); ungetc(c,inFILE);
  if(c=='0') return -1; if((c<'0') || ('9'<c)) return 0; else return 1;
}

void Print_PPL(PolyPointList *_P, const char *comment){
  int i,j;
  if(_P->np>20){
    fprintf(outFILE,"%d %d  %s\n",_P->np,_P->n,comment);
    for(i=0;i<_P->np;i++) {
      for(j=0;j<_P->n;j++) fprintf(outFILE,"%d ",(int) _P->x[i][j]); 
      fprintf(outFILE,"\n");}}
  else {
    fprintf(outFILE,"%d %d  %s\n",_P->n,_P->np,comment);
    for(i=0;i<_P->n;i++) {
      for(j=0;j<_P->np;j++) fprintf(outFILE," %4d",(int) _P->x[j][i]); 
      fprintf(outFILE,"\n");}}
}

void Print_VL(PolyPointList *_P, VertexNumList *_V, const char *comment){
  int i,j;
  if(_V->nv>20){
    fprintf(outFILE,"%d %d  %s\n",_V->nv,_P->n,comment);
    for(i=0;i<_V->nv;i++) {
      for(j=0;j<_P->n;j++) fprintf(outFILE,"%d ",(int) _P->x[_V->v[i]][j]); 
      fprintf(outFILE,"\n");}}
  else {
    fprintf(outFILE,"%d %d  %s\n",_P->n,_V->nv,comment);
    for(i=0;i<_P->n;i++) {
      for(j=0;j<_V->nv;j++) fprintf(outFILE," %4d",(int) _P->x[_V->v[j]][i]); 
      fprintf(outFILE,"\n");}}
}

void Print_EL(EqList *_E, int *n, int suppress_c, const char *comment){
  int i,j;
  fprintf(outFILE,"%d %d  %s\n",_E->ne,*n,comment);
  for(i=0;i<_E->ne;i++) {
    for(j=0;j<*n;j++) fprintf(outFILE," %3d",(int) _E->e[i].a[j]); 
    if (!suppress_c) fprintf(outFILE," %5d",(int) _E->e[i].c);
    fprintf(outFILE,"\n");}
}

void Print_Matrix(Long Matrix[][VERT_Nmax], int n_lines, int n_columns, 
              const char *comment){
  int i,j;
  fprintf(outFILE,"%d %d  %s\n",n_lines, n_columns, comment);
  for(i=0;i<n_lines;i++) {
    for(j=0;j<n_columns;j++) fprintf(outFILE," %3d",(int) Matrix[i][j]); 
    fprintf(outFILE,"\n");}
}
int  auxString2Int(char *c,int *n)
{    int j=0; *n=0;  while(('0'<=c[j])&&(c[j]<='9')) *n = 10*(*n)+c[j++]-'0';
     if(j) while(c[j]==' ') j++; return j;
}
void CWSZerror(char *c)
{    printf("Format error %s in Read_CWS_Zinfo\n",c);exit(0);
}
void Print_CWS_Zinfo(CWS *CW)
{    int i,j; if(CW->nw) for(i=0;i<CW->nz;i++) 
     {      fprintf(outFILE,"/Z%d: ",CW->m[i]);
      for(j=0;j<CW->N;j++) fprintf(outFILE,"%d ",CW->z[i][j]);
     }      
}
int  Read_CWS_Zinfo(FILE *inFILE,CWS *CW)             /* return !EOF */
{    int *nz=&CW->nz; int i=0,n; char c[999],b=' '; *nz=0; 
     for(n=0;n<999;n++)
     {      c[n]=fgetc(inFILE); if(feof(inFILE)) return 0; if(c[n]=='\n') break;
     }  if(n==999) {puts("Out of space in Read_CWS_Zinfo");exit(0);}
     while(c[i]==b)i++; if((c[i]=='=')&&(c[i+1]=='d'))i+=2; while(c[i]==b)i++; 
     /*     printf("Read_CWS_Zinfo n=%d -> ",n);
      {int I; for(I=0;I<n;I++)printf("%c",c[I]);puts("");} */
     while(i<n)
     {      int j, k, s; if((c[i]!='/')||(c[i+1]!='Z')) return 1;
        i+=2; assert(*nz<POLY_Dmax);
        if((j=auxString2Int(&c[i],&CW->m[*nz])))
      {   if(c[i+j]!=':') CWSZerror(":"); else {i+=j+1; while(c[i]==b)i++;}
          for(k=0;k<CW->N;k++)
          { if((j=auxString2Int(&c[i],&CW->z[*nz][k]))) 
                  {if((j)) i+=j; else CWSZerror("missing");} /* ????? */
          }   s=0;
          for(k=0;k<CW->N;k++) s+=CW->z[*nz][k]; 
          /*      if(s % CW->m[*nz])  CWSZerror("det!=1"); */
          (CW->nz)++;
        }   else CWSZerror("Order");
     }      return 1;
}

int  Read_CWS_PP(CWS *_CW, PolyPointList *_P)
{    int i, j, FilterFlag=(inFILE==NULL); 
     int IN[AMBI_Dmax*(AMBI_Dmax+1)], S; 
     static int InputOK;
     _CW->nw=_CW->N=_CW->nz=0;
     if(FilterFlag) inFILE=stdin; 
     else if(inFILE==stdin)     {   
       printf("Degrees and weights  `d1 w11 w12 ... d2 w21 w22 ...'\n  or `#");
       printf("lines #colums' (= `PolyDim #Points' or `#Points PolyDim'):\n");}
     for(i=0;i<AMBI_Dmax*(AMBI_Dmax+1);i++)
     { char c;
       while(' '==(c=fgetc(inFILE ))); ungetc(c,inFILE);    /* read blanks */
       if(IsNextDigit()) fscanf(inFILE,"%d",&IN[i]); else break;
     }
     if(i==0) { if(!InputOK){puts("-h gives you help\n"); exit(0);}
      else return 0; } InputOK++;
     if(i==1) { puts("Error in INPUT: need at least 2 numbers!"); exit(0);}
     if(i==2)                                   /* READ PolyPointList */
     {      int tr=0;    while('\n'!=fgetc(inFILE));      /* read to end of line */
        if(IN[0]>IN[1]) {tr=IN[0];IN[0]=IN[1];IN[1]=tr;} tr=!tr;
      if(IN[0]>POLY_Dmax) {puts("increase POLY_Dmax!"); exit(0);}
      if(IN[1]>POINT_Nmax) {puts("increase POINT_Nmax!"); exit(0);}
      if((inFILE==stdin)&&!FilterFlag)  
        printf("Type the %d coordinates as %s=%d lines with %s=%d colums:\n",
            IN[0]*IN[1], tr ? "dim" : "#pts", tr ? IN[0] : IN[1],
            tr ? "#pts" : "dim", tr ? IN[1] : IN[0]);
      _P->n=IN[0]; _P->np=IN[1];
/* allow all numbers in one string (distributed over lines) or
 * matrix blocks with trailing comments   */
      if(tr) for(i=0;i<IN[0];i++) for(j=0;j<IN[1];j++)
        {   int X; fscanf(inFILE,"%d",&X); _P->x[j][i]=X;   }
      else for(i=0;i<IN[1];i++) for(j=0;j<IN[0];j++)
        {   int X; fscanf(inFILE,"%d",&X); _P->x[i][j]=X;   }
      /* Finish_Poly_Points(_P); */
      while(fgetc(inFILE )-'\n') if(feof(inFILE)) return 0;/* read to EOL */
      if(FilterFlag) inFILE=NULL; return 1;
     }  assert(i!=3);
     S=IN[i-1]; for(j=0;j<i-1;j++) if(S<IN[j]) break; 
     if(j==i-1)                      /* Single Weights with: "w1 ... d" */
     {      _CW->nw=1; _CW->d[0]=S; _CW->N=i-1; 
        assert(_CW->N <= POLY_Dmax+1);                /* Increase POLY_Dmax */ 
      for(j=0;j<_CW->N;j++) {_CW->W[0][j]=IN[j]; assert(IN[j]>0);} 
      goto MAP;
     }
     S=_CW->d[_CW->nw]=IN[0];                             /* ASSIGN WEIGHTS */
     for(j=1; j<i; j++)
     {      if(S<IN[j]) 
      {   if(S) { puts("Error in INPUT: degree vs. weights!"); exit(0);}
          else  
          { S=_CW->d[++(_CW->nw)]=IN[j];
            if(1==_CW->nw) 
            {   if(i%(1+_CW->N)) {puts("INPUT error: numbers?"); exit(0);}
            }   if(j%(1+_CW->N)) {puts("INPUT error: degrees?"); exit(0);}
          }
      }
      else
      {   _CW->W[_CW->nw][j - 1 - (1 + _CW->N)*_CW->nw] = IN[j]; S-=IN[j];
          if(0==_CW->nw) if(++(_CW->N)>AMBI_Dmax) 
          {   puts("Increase AMBI_Dmax!"); exit(0);
          }
      }
     }
     ++(_CW->nw);
     if(_CW->N-_CW->nw > POLY_Dmax){
       printf("Please increase POLY_Dmax to at least %d\n", _CW->N-_CW->nw);
       exit(0);} /* increase POLY_Dmax */

MAP: for(i=0;i<_CW->nw;i++)               /* check consistency of CWS */
     {      Long sum=_CW->d[i], *w=_CW->W[i];
      for(j=0;j<_CW->N;j++) { assert(w[j]>=0); sum-=w[j]; }
      if(sum){ printf("Use poly.x -w for (single) WeightSystems with ");
            printf("d!=\\sum(w)\n(only Read_Weight makes the correct ");
            puts("PolyPointList in that case"); exit(0);}
     }
     _CW->nz=0;
     if(!Read_CWS_Zinfo(inFILE,_CW))                     /* read Z to EOL */
     {      if(!InputOK) puts("-h gives you help\n"); return 0;}
     Make_CWS_Points(_CW, _P);                           /* now make POLY */
     if(FilterFlag) inFILE=NULL; return 1;
}

void Print_CWH(CWS *_W, BaHo *_BH){
  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);
  if(_BH->np)     {
    fprintf(outFILE,"M:%d %d N:%d %d",_BH->mp,_BH->mv,_BH->np,_BH->nv);
    if (_BH->n == 3) fprintf(outFILE," Pic:%d Cor:%d",_BH->pic, _BH->cor); 
    if (_BH->n > 3) {
      fprintf(outFILE," H:%d",_BH->h1[1]); 
      for(i=2;i<_BH->n-1;i++) fprintf(outFILE,",%d",_BH->h1[i]);
      if(_BH->n==4) fprintf(outFILE," [%d]",2*(_BH->h1[1]-_BH->h1[2])); 
      if(_BH->n==5) 
      fprintf(outFILE," [%d]",48+6*(_BH->h1[1]-_BH->h1[2]+_BH->h1[3]));}}
  else if(_BH->mp)
    fprintf(outFILE,"M:%d %d F:%d",_BH->mp,_BH->mv,_BH->nv);
  else fprintf(outFILE,"V:%d F:%d",_BH->mv,_BH->nv);
  fprintf(outFILE,"\n");
}

/*  ==========          END of I/O functions          ==========  */


/*  ==========       Solve W-Eq. -> triangular LatticeBasis
 *   
 *   B[0]:= (-n1, n0, 0, ...) / GCD(n0,n1);
 *   B[i]:= (0,...,0,g/G,0,...)- (ni/G) * ExtGCD.K(n0,...,n(i-1),0,...);
 *        with g=GCD(n0,...,n[i-1]); G=GCD(g,ni);
 *
 *   CWS by iteration: NextWeight=W[]*B; NewB=Basis(NextW); B[i+1]=B[i]*NewB;
 *                                                        */

void PrintBasis(CWLatticeBasis *_B)
{    int i,j;                puts("Basis:");
     for(i=0;i<_B->n;i++)
     {      for(j=0;j<_B->N;j++) fprintf(outFILE,"%6d ",(int) _B->x[i][j]); 
      puts("");
     }      puts("End of Basis  - -");
}

void Orig_Solve_Next_WEq(Long *NW, CWLatticeBasis *_B)
{    int i, j, P=0, p[AMBI_Dmax]; Long W[AMBI_Dmax], G;     _B->n=_B->N-1; 
     for(i=0;i<_B->N;i++) 
     {      for(j=0;j<_B->n;j++) _B->x[j][i]=0;                   /* init B.x=0 */
      if(NW[i]) {p[P]=i; W[P++]=NW[i];}         /* non-zero weights */
     }
     if(P<2) puts("need two non-zero weights in  >>Solve_Next_WEq<<");
     for(i=0;i<p[0];i++) _B->x[i][i]=1; while((++i)<p[1]) _B->x[i-1][i]=1;
     G=Fgcd(W[0],W[1]); if(W[0]/G<0) G=-G;
     _B->x[i-1][p[0]]=-W[1]/G; _B->x[i-1][p[1]]=W[0]/G;
     j=2; while(++i<_B->N)
     {      if(NW[i])
      {   int k; Long *X=_B->x[i-1], K[AMBI_Dmax], g=REgcd(W, &j, K); 
          G=Fgcd(g,NW[i]); if(g/G<0) G=-G; X[i]= g/G; g=W[j]/G;
          for(k=0;k<j;k++) X[p[k]]=-K[k]*g; j++;
      }
      else _B->x[i-1][i]=1;
     }
}
void Solve_Next_WEq(Long *NW, CWLatticeBasis *_B)
{    Long W[AMBI_Dmax], *X[AMBI_Dmax], GLZ[AMBI_Dmax][AMBI_Dmax];
     int i, j, P=0, p[AMBI_Dmax]; _B->n=_B->N-1;
#ifdef TEST_Wbase
     Orig_Solve_Next_WEq(NW,_B); PrintBasis(_B);
#endif
     for(i=0;i<_B->N;i++) 
     {      for(j=0;j<_B->n;j++) _B->x[j][i]=0;             /* init B.x=0 */
      if(NW[i]) {p[P]=i; X[P]=GLZ[P]; W[P++]=NW[i];}  /* non-zero weights */
     }      if(P>1) W_to_GLZ(W,&P,X);           /* P>1, compute GLZ */
     else {/* printf("P=%d W[0]=%d for W_to_GLZ\n",P,W[0]);exit(0);*/assert(P);
      for(i=0;i<p[0];i++)_B->x[i][i]=1; while((++i)<_B->N)_B->x[i-1][i]=1;
      return; }
     for(i=1;i<P;i++) if(X[i][i]<0) for(j=0;j<=i;j++) X[i][j] *= -1;
     for(i=0;i<p[0];i++) _B->x[i][i]=1; while((++i)<p[1]) _B->x[i-1][i]=1;
     _B->x[i-1][p[0]]=X[1][0]; _B->x[i-1][p[1]]=X[1][1]; j=2;
     while(++i<_B->N)
     {      if(NW[i])
      {   int k; Long *B=_B->x[i-1]; for(k=0;k<=j;k++) B[p[k]]=X[j][k]; j++;
      }
      else _B->x[i-1][i]=1;
     }
#ifdef TEST_Wbase
      /* for(i=0;i<_B->N;i++)printf("%d ",NW[i]);puts("=NW");for(i=0;i<P;i++)
      {for(j=0;j<P;j++)printf("%4d ",GLZ[i][j]);puts("=GL");} */
      printf("New version: "); PrintBasis(_B);
#endif
#ifdef USE_Old_Wbase
      Orig_Solve_Next_WEq(NW,_B);
#endif
}
void Make_CWS_Basis(CWS *_C, CWLatticeBasis *_B){    
  int i,j,k,l; Long W[AMBI_Dmax]; CWLatticeBasis NB, AuxB; 
  AuxB.N=_B->N=_C->N;
  Solve_Next_WEq(_C->W[0],_B);
  for(i=1;i<_C->nw;i++){ 
    for(j=0;j<_B->n;j++){   
      W[j]=0; 
      for(k=0;k<_B->N;k++) W[j]+=_C->W[i][k]*_B->x[j][k]; }
    NB.N=_B->n; 
    Solve_Next_WEq(W,&NB);
    AuxB.n=NB.n;
    for(j=0;j<AuxB.N;j++) for(k=0;k<AuxB.n;k++) {   
      AuxB.x[k][j]=0; 
      for(l=0;l<NB.N;l++) AuxB.x[k][j] += NB.x[k][l] * _B->x[l][j]; }
    *_B=AuxB;     }
  /* init coordinate hyperplanes */
  _C->B.ne=_B->N;   /* == transposed of Basis */
  for(i=0;i<_B->N;i++) {
    for(j=0;j<_B->n;j++) _C->B.e[i].a[j]=_B->x[j][i];
    _C->B.e[i].c=1;}
}

void Poly_To_Ambi(CWLatticeBasis *_B, Long *x, Long *X)
{    int i, j;
     for(i=0;i<_B->N;i++)
     {      X[i]=1; for(j=0;j<_B->n;j++) X[i] += _B->x[j][i] * x[j];
     }
}
/*   B.x[i][A] is lower triangular: zero for i<p and Amin[p] <= A <Amin[p+1]
 *   X^A=x^p B_p^A+1, hence  for A \in [ Amin[p],Amin[p+1] ) and all k:
 *   x^iB_i^A \in [0,d_k/w_k^A] - 1 - \sum_{l<i} x^l B_l^A
 *                                                         */

void QuotZ_2_SublatG(Long Z[][VERT_Nmax],int *zm,Long *M,int *d,
     Long G[][POLY_Dmax]);                /* normalize and diagonalize Z */
void CWS_2_SublatZ(CWS *C,CWLatticeBasis *B,                       /* in */
      int *m, Long *M, Long G[POLY_Dmax][POLY_Dmax])              /* out */
{    int i,j,k,d=C->N-C->nw; Long Z[POLY_Dmax][VERT_Nmax]; *m=C->nz; 
     for(i=0;i<C->nz;i++){ M[i]=C->m[i]; for(j=0;j<d;j++)
     {      Z[i][j]=0; for(k=0;k<C->N;k++) Z[i][j]+=C->z[i][k]*B->x[j][k];}
     }      /* PrintBasis(B); */    QuotZ_2_SublatG(Z,m,M,&d,G);
}
void Reduce_PPL_2_Sublat(PolyPointList *P,int *nm,Long *M,Long G[][POLY_Dmax])
{    int i,j,n=0,N; for(N=0;N<P->np;N++) 
     {      Long X[POLY_Dmax]; for(i=0;i<P->n;i++)
      {   X[i]=0; for(j=0;j<P->n;j++) X[i]+=G[i][j]*P->x[N][j];
      }   for(i=0;i<*nm;i++) if(X[i]%M[i]) break; if(i<*nm) continue;
      for(i=0;i<*nm;i++) P->x[n][i]=X[i]/M[i]; 
      for(i=*nm;i<P->n;i++) P->x[n][i]=X[i]; n++;
     }      P->np=n; /* Print_PPL(P,"");printf("nm=%d\n",*nm); */
}
Long PD_Floor(Long N,Long D)  /*  assuming PosDenom  D>0:  F <= N/D < F+1  */
{    Long F=N/D; return (F*D>N) ? F-1 : F; 
}
void Make_CWS_Points(CWS *Cin, PolyPointList *_P)
{    int i, j, Amin[POLY_Dmax+1]; Long *x=_P->x[_P->np=0], xmin[POLY_Dmax], 
      xmax[POLY_Dmax], Xmax[AMBI_Dmax], xaux[POLY_Dmax], L, R; CWS *_C=Cin;
     CWLatticeBasis B; Long G[POLY_Dmax][POLY_Dmax],M[POLY_Dmax];int m=Cin->nz;
#ifndef NO_COORD_IMPROVEMENT        /* ==== Perm Coord Improvement ==== */
     int pi[AMBI_Dmax]; CWS Caux; _C=&Caux; CWS_to_PermCWS(Cin,_C, pi);
#endif                               /* = End of Perm Coord Improvement = */
     Make_CWS_Basis(_C, &B);               /* make `triangular' Basis */
#ifndef NO_COORD_IMPROVEMENT        /* ==== Perm Coord Improvement ==== */
     assert(_C->N==_C->B.ne);
     for(i=0;i<Cin->N;i++) Cin->B.e[i]=_C->B.e[pi[i]]; Cin->B.ne=_C->N;
#endif                               /* = End of Perm Coord Improvement = */
     i=_P->n=B.n; Amin[0]=0; Amin[B.n]=j=B.N;         /* inversion structure */
     while(--i) {while(!B.x[i-1][--j]); Amin[i]=++j;}
     for(i=0;i<B.N;i++)                              /* compute Xmax */
     {      Xmax[i]=0; for(j=0;j<_C->nw;j++) if(_C->W[j][i])
      {   L=_C->d[j]/_C->W[j][i]; 
          if(Xmax[i]) {if(L<Xmax[i]) Xmax[i]=L;}
          else Xmax[i]=L;
      }
     }
     j=B.n-1; i=Amin[j+1]-1; R=B.x[j][i];     /* compute xmin[j] and xmax[j] */
     xmin[j]=-PD_Floor(1,R); xmax[j]=PD_Floor(Xmax[i]-1,R);   /* since R > 0 */
/** /printf("R=%2d:  %2d <= x[%d=&%d] <= %2d\n",R,xmin[j],j,i,xmax[j]);/ **/
     while((i--)>Amin[j])     /* if(R=B.x[B.n-1][i]):  R!=0  => new limits */
     {      Long Low=-1, Upp=Low+Xmax[i];               R=B.x[B.n-1][i];
      if(R>0) {if(xmax[j]>(L=PD_Floor(Upp,R)))   xmax[j]=L;
             if(xmin[j]<(L=-PD_Floor(-Low,R))) xmin[j]=L;}
      else  {if(xmax[j]>(L=PD_Floor(-Low,-R))) xmax[j]=L;
             if(xmin[j]<(L=-PD_Floor(Upp,-R))) xmin[j]=L;}
/** /printf("R=%2d:  %2d <= x[%d=&%d] <= %2d\n",R,xmin[j],j,i,xmax[j]);/ **/
     }                           /* this completes the limits for x[B.n-1] */
     x[j]=xmin[j]; while(j<B.n)
     {      int k; if(x[j]>xmax[j]) {if(B.n==(++j)) break; else x[j]++;}
      else        /* compute limits[j-1] and initialize x[j-1] */
      {   Long Low=-1, Upp=Xmax[i=Amin[j--]-1]; int RangeFlag=0;
          for(k=j+1;k<B.n;k++) Low-=x[k]*B.x[k][i];      /* compute offset */
          Upp+=Low; R=B.x[j][i]; 
          xmin[j]=-PD_Floor(-Low,R); xmax[j]=PD_Floor(Upp,R);
/** /printf("R=%2d:  %2d <= x[%d=&%d] <= %2d\n",R,xmin[j],j,i,xmax[j]);/ **/
          while((i--)>Amin[j]) if((R=B.x[j][i]))    /* R!=0  => new limits */
          { Low=-1; Upp=Xmax[i];
            for(k=j+1;k<B.n;k++) Low-=x[k]*B.x[k][i];
            Upp+=Low; R=B.x[j][i]; 
            if(R>0) {if(xmax[j]>(L=PD_Floor(Upp,R)))   xmax[j]=L;
                   if(xmin[j]<(L=-PD_Floor(-Low,R))) xmin[j]=L;}
            else  {if(xmax[j]>(L=PD_Floor(-Low,-R))) xmax[j]=L;
                   if(xmin[j]<(L=-PD_Floor(Upp,-R))) xmin[j]=L;}
/** /printf("R=%2d:  %2d <= x[%d=&%d] <= %2d\n",R,xmin[j],j,i,xmax[j]);/ **/
          }                   /* completes limits for x[] except */
          else                /*   when R=0 and X out of Range:  */
          { Long X=1; for(k=j+1;k<B.n;k++) X+=x[k]*B.x[k][i];
            if((X<0)||(X>Xmax[i])) RangeFlag=1;
          }
          if(RangeFlag) ++x[++j]; else x[j]=xmin[j];            
          if(j==0)
          {   while(x[0]<=xmax[0]) 
            {   Long *y; 
                if((++_P->np) < POINT_Nmax) 
                { y=(_P->x[_P->np]); for(k=0;k<B.n;k++) y[k]=x[k]; 
                }
                else if(_P->np == POINT_Nmax)
                { y=xaux; for(k=0;k<B.n;k++) y[k]=x[k];
                }
                else {puts("Increase POINT_Nmax");exit(0);}
                x=y; ++x[0];
            }
              x[j=1]++;
          }
      }
     }      if(m)CWS_2_SublatZ(_C,&B,&m,M,G);if(m)Reduce_PPL_2_Sublat(_P,&m,M,G);
}

/*  ==========    Coordinate improvement via CWS-Permutations   ==========  */

typedef struct {Long s, g, *N, **G; int *p;}          Tri_GLZ_MPaux;  
Long Tri_GLZ_Norm(int *d, Long **S)             /* sum-norm = default */
{    Long norm=0; int i,j; for(i=0;i<*d;i++)for(j=0;j<*d;j++)
#ifdef      MAX_NORM                                  /* max(|*|) */
     {      Long x=S[i][j]; if(x<0) x=-x; if(x>norm) norm=x; } return norm;
#else                                           /* sum(|*|) */
     {      Long x=S[i][j]; norm += (x>0) ? x : -x ; } return norm; 
#endif
}
void Tri_GLZ_Basis_Perm(int *d,int *pi, /* int *pinv, */ Tri_GLZ_MPaux * AP)
{    Long g, norm, N[AMBI_Dmax], M[AMBI_Dmax][AMBI_Dmax], *S[AMBI_Dmax];  
     int i, j; /* if(pi[0]>pi[1]) return; *//* eliminate equiv. permut. */
     for(i=0;i<*d;i++) { N[i]=AP->N[pi[i]]; S[i]=M[i];}
     g=W_to_GLZ(N,d,S); norm=Tri_GLZ_Norm(d,S);
     if(AP->s) assert(g==AP->g); else AP->g=g;
#ifdef      TEST_MIN_GLZ
     {      Long err=0, max=0, tsum=0, pos; for(i=0;i<*d;i++)
      {   for(j=0;j<*d;j++) {pos=S[i][j]; if(pos<0)pos=-pos; tsum+=pos; 
          if(max<pos)max=pos;} if(err==0) { 
            for(j=0;j<*d;j++) err+=N[j]*S[i][j]; if(i==0)err-=g;}
      }   assert((norm==tsum)||(norm==max)); 
      printf("max=%lld, sum=%lld\n",(long long) max, (long long) tsum);
      if(err) { for(i=0;i<*d;i++) printf("%d",pi[i]); 
      printf("  g=%d  norm=%d",(int)g,(int)norm); puts(""); 
      for(i=0;i<*d;i++){printf("S[%d]=",i);for(j=0;j<*d;j++)printf(" %5d",
      S[i][j]);if(!i)printf("   g=%d  norm=%d",(int)g,(int)norm);puts("");}}
     }
#endif
     if((0 == AP->s) || (norm < AP->s))            /* init or improve AP->G */
     {      for(i=0;i<*d;i++)for(j=0;j<*d;j++)AP->G[i][pi[j]]=S[i][j]; AP->s=norm;
      for(i=0;i<*d;i++) AP->p[i]=pi[i];
     }
}

#ifndef NO_COORD_IMPROVEMENT  
                        /* improved by permutation pi=pinv^-1 */
Long Wperm_to_GLZ(Long *W, int *d, Long **G, int *P)
{    Tri_GLZ_MPaux AS; int i, j, pi[AMBI_Dmax], pinv[AMBI_Dmax];
     AS.s=0; AS.N=W; AS.G=G; AS.p=P;
     Map_Permut(d,pi, /* pinv,*/ (Tri_GLZ_Basis_Perm),(void *) &AS);
     for(i=0;i<*d;i++)for(j=0;j<*d;j++)G[i][j]=AS.G[i][j];
#ifdef      TEST_MIN_GLZ
      for(i=0;i<*d;i++){printf("  G[%d]=",i);for(j=0;j<*d;j++)printf(" %5d",
      G[i][j]);if(!i)printf("  g=%d norm=%d",(int)AS.g,(int)AS.s);puts("");}
#endif
     return AS.g;
}
void CWS_to_PermCWS(CWS *Cin, CWS *C, int *pi)
{    int i, j, N=C->N=Cin->N, n=C->nw=Cin->nw, l=0, A[AMBI_Dmax]; 
     Long *X, *G[AMBI_Dmax], M[AMBI_Dmax][AMBI_Dmax], W[AMBI_Dmax]; 
     for(j=1;j<n;j++) if(Cin->d[j]>Cin->d[l]) l=j; X=Cin->W[l]; 
     j=0; for(i=0;i<N;i++) if(X[i]) {A[j]=i; W[j]=labs(X[i]); G[j]=M[j]; j++;}
     Wperm_to_GLZ(W,&j,G,pi); for(i=0;i<j;i++) pi[i]=A[pi[i]];
     for(i=0;i<N;i++) if(X[i]==0) pi[j++]=i;
     C->d[0]=Cin->d[l]; for(i=0;i<N;i++) C->W[0][i]=Cin->W[l][pi[i]];
     for(j=1;j<n;j++) {int L=j-(j<=l); C->d[j]=Cin->d[L];
      for(i=0;i<N;i++) C->W[j][i]=Cin->W[L][pi[i]];}
     C->nz=Cin->nz; for(j=0;j<C->nz;j++){C->m[j]=Cin->m[j];
      for(i=0;i<N;i++) C->z[j][i]=Cin->z[j][pi[i]];}
}
#endif

Generated by  Doxygen 1.6.0   Back to index