pelclyal.cc

Gambit 是一个游戏库理论软件

/* cly_all.c */

#include "pelclyal.h"


/************************************************************************/
/****************** implementation code from cly_cells.c ****************/
/************************************************************************/

/*
******************************************************************
** Global  Variables
******************************************************************
*/
int cly_Npts=0;      /* the number of points in the config */
int cly_N=0;        /* the dimension of the cayley point config */
int cly_R=0;        /* the number of point configs */
int cly_Dim=0;        /* the dimension of the Aset */
int next_id=1; /*for debuging cells are labled by a unique id*/

/*
** The matrices here are used for now to avoid problems,
** they will eventually be stored with the individual cells.
*/

HMatrix cly_U=0;   /* factor matrix */ 
HMatrix cly_M=0;   /* an n+1xn+1 matrix   */
HMatrix cly_L=0;   /* an n+1 vector   */
Imatrix cly_T=0;   /* an R vector */

/*
** two temporary variables used by the Hint macroes
*/
Hint cly_temp;
int cly_det;
/*
** controll parameters
**
*/
int cly_order=TRUE;
int cly_lift=TRUE;


/* 
** Display functions: mainly for debugging
**  Ipnt_print  print complete list of data fields.
**  Ipnt_list_print  print list print lables of all points accesable
**                   from a given points through nex pointers
*/
void Ipnt_fprint(FILE *fout,Ipnt p){
  int i;
#ifdef LOG_PRINT
  fprintf(fout,"%% %s,",Ipnt_lable(p))
#endif
;
  for (i=1;i<=cly_N;i++) 
#ifdef LOG_PRINT
fprintf(fout,"%d, ",Ipnt_coord(p,i))
#endif
;
#ifdef LOG_PRINT
  fprintf(fout,"lift = %d, idx =%d \n",Ipnt_lift(p),Ipnt_idx(p))
#endif
;
 }

void Ipnt_list_fprint(FILE *fout,Ipnt p){
#ifdef LOG_PRINT
 fprintf(fout,"%%{Pl\n")
#endif
;
 for(;p!=0;p=Ipnt_next(p)) Ipnt_fprint(fout,p);
#ifdef LOG_PRINT
 fprintf(fout,"%%Pl}\n")
#endif
;
}

/*
** Creator/Destructor
*/
 Ipnt Ipnt_new(node pnt,int r){
    Ipnt p;
    int j;
    p=(Ipnt)mem_malloc(sizeof(struct Ipnt_t));
    if (p==0) bad_error("malloc 1 failed in new_Ipnt");
    p->point=pnt;
    p->idx=r;
    p->lable=(char *)mem_strdup(pnt_lable(pnt));
    p->coords=(int *)mem_malloc(cly_N*sizeof(int));
    if (p==0) bad_error("malloc 2 failed in new_Ipnt");
    for(j=1;j<=cly_Dim;j++) Ipnt_coord(p,j)=*IVref(pnt_coords(pnt),j);
    for(j=cly_Dim+1;j<=cly_N;j++) Ipnt_coord(p,j)=0;
    if (r>1) Ipnt_coord(p,cly_Dim+r-1)=1;
    p->lift=0;
    p->next=0;
    return p;
}


void Ipnt_free(Ipnt p){
  if (p!=0){
     mem_free(p->lable);
     mem_free((char*)(p->coords));
     mem_free((char *)p);
  }
 }

void points_free(Ipnt p){
  Ipnt jk;
  while(p!=0){
    jk=p;
    p=Ipnt_next(p);
    Ipnt_free(jk);
  }
 }

void lift_original_points(Ipnt p){
 for( ;p!=0; p=Ipnt_next(p))
     *IVref(pnt_coords(Ipnt_pnt(p)),cly_Dim+1)=Ipnt_lift(p);
}

/**************************************************************
** functions for basic acces and manipulation of cells.
***************************************************************/

/*
** cell_new(int n,int r)
**   reserve space for a new cell, for holding an n-simplex
**   the points and pointers fields are initialized to null
*/
 cell cell_new(int n,int r){
  cell c;
  int i;
  if ((c=(cell)mem_malloc(sizeof(struct cell_t)))==0)
                bad_error("malloc failed in cell_new()");
  if ((c->points=(Ipnt *)mem_malloc((n+1)*sizeof(Ipnt)))==0)
                bad_error("malloc failed in cell_new()");
  if ((c->ptrs=(cell *)mem_malloc((n+1)*sizeof(cell)))==0)
                bad_error("malloc failed in cell_new()");
  c->H=HMnew(n,n+1);
  c->U=HMnew(n,n);
  c->T=Ivector_new(r);
  c->norm=HVnew(n+1);
  c->n=n;
  c->r=r;
  c->id=next_id++;
  c->freeptrs=n+1;
  c->volume=0;
  c->Tindx=0;
  c->next=0;
  for(i=0;i<=n;i++) {c->ptrs[i]=0; c->points[i]=0;}
  return c;
}

/*
** cell_free(cell c)
**           free space allocated for a cell
**           (does not effect any of the objects pointed to by the
**            points or pointers fields).
*/
 void cell_free(cell c){
  if (c!=0){
     mem_free((char *)c->points);
     mem_free((char *)c->ptrs);
     HMfree(c->norm);
     Imatrix_free(c->T);
     HMfree(c->U);
     HMfree(c->H);
     mem_free((char *)c);
  }
}
 void subdiv_free(cell c){
 cell jk;
 while(c!=0){ jk=c; c=cell_next(c); cell_free(jk);}
}

/*
**   cell_print(cell c)-   print the contents of a single cell.
*/
 void cell_fprint(FILE *fout, cell c){
   int i,k,tog;

#ifdef LOG_PRINT
   fprintf(fout,"%% {")
#endif
;
   for(k=1;k<=cly_R;k++){
#ifdef LOG_PRINT
     fprintf(fout,"{")
#endif
;
     tog=0;
     for (i=0;i<=c->n;i++){
       if (Ipnt_idx(cell_pnt(c,i))==k){
         if (tog==0) tog=1;
         else printf(",");
     #ifdef LOG_PRINT
    fprintf(fout," %s",Ipnt_lable(cell_pnt(c,i)))
#endif
;
       }
     }
#ifdef LOG_PRINT
     fprintf(fout,"}")
#endif
;
     if (k<cly_R) fprintf(fout,",");
   }
#ifdef LOG_PRINT
   fprintf(fout,"}   < %d",cell_type(c,1))
#endif
;
   for(i=2;i<=cly_R;i++) 
#ifdef LOG_PRINT
fprintf(fout,", %d",cell_type(c,i))
#endif
;
#ifdef LOG_PRINT
   fprintf(fout," >        vol=%d",cell_volume(c));
   fprintf(fout,"\n")
#endif
;
}

/*
** subdiv_print(cell c)-   print contents of all cells in a
**    subdivision, a list of cells linked through their next fields
*/
 void subdiv_fprint(FILE *fout, cell c){
#ifdef LOG_PRINT 
  fprintf(fout,"\n")
#endif
;
  while(c!=0){
      cell_fprint(fout,c);
      c=cell_next(c);
  }

}

/*
** point_is_in(Ipnt pt, cell c)
**   Input:  a point pt,
**           a cell c.
**   Output: TRUE if cell contains pt
**           FALSE otherwise.
** (a point is determined by its memory address, points with same
**  coordinates and lables but stored in diferent places are still
**  distinct).
*/
 int point_is_in(Ipnt pt, cell c){
  int i;
  for(i=0;i<=cell_n(c);i++)
          if (pt==cell_pnt(c,i)) return TRUE;
  return FALSE;
 }
 int cell_type_cmp(cell c1,cell c2){
  int i,t;
  if (c1==0 && c2==0) return 0;
  else if (c1==0) return 1;
  else if (c2==0) return -1;
  for(i=1;i<=cly_R;i++){
     t=cell_type(c1,i)-cell_type(c2,i);
     if (t!=0) return t;
  }
  return 0;
}

 void order_subdiv(cell *Subdiv){
  cell Res=0,S=*Subdiv, ptr=0,pts=0;
  while(S!=0){
    pts=S;
    S=cell_next(S);
    if (Res==0 || cell_type_cmp(pts,Res)>=0){
      cell_next(pts)=Res;
      Res=pts;
    }
    else {
      ptr=Res;
      while(cell_next(ptr)!=0 && cell_type_cmp(pts,cell_next(ptr))<0){
         ptr=cell_next(ptr);
      }
      cell_next(pts)=cell_next(ptr);
      cell_next(ptr)=pts;
    }
  }
 *Subdiv=Res;
}
 int fprint_all_volumes(FILE *fout,cell S){
  int vol=0, i;
  while(S!=0){
     vol+=cell_volume(S);
     if (cell_next(S)==0 || cell_type_cmp(S,cell_next(S))!=0){
      #ifdef LOG_PRINT
  fprintf(fout,"%% Vol( %d", cell_type(S,1))
#endif
;
        for(i=2;i<=cly_R;i++) 
#ifdef LOG_PRINT
fprintf(fout,", %d",cell_type(S,i))
#endif
;
#ifdef LOG_PRINT
        fprintf(fout,")= %d\n",vol)
#endif
;
        vol=0;
     }
     S=cell_next(S);
  }
  return 1;
}

 int cell_set_volume(cell c){
 int i,j,r,row;

 cly_T=Imatrix_resize(cly_T,1,cly_R);
 cly_U=HMresize(cly_U,cly_Dim,cly_Dim);
 cly_M=HMresize(cly_M,cly_Dim,cly_Dim);

 for(i=1;i<=cly_R;i++)*IVref(cly_T,i)=-1;

 row=1;
 for(i=0;i<=cly_N;i++){
   r=Ipnt_idx(cell_pnt(c,i));
   if (*IVref(cly_T,r)==-1){
             *IVref(cly_T,r)=i;
   }
   else{
    for(j=1;j<=cly_Dim;j++){
      HLMset(cly_M,row,j,cell_point(c,i,j)-cell_point(c,*IVref(cly_T,r),j));
    }
    row++;
   }
 }
 HMdet(cly_M,cly_U,r);
 return (cell_volume(c)=r);
}

/*
** subdiv_union
**    Input:        subdivisions SD1 and SD2
**   Output:        SD2 appended to SD1.
**   Side Effects:  if SD1 is not null SD1 also points to union.
*/
 cell cly_subdiv_union(cell SD1, cell SD2){
  cell ptr=SD1;
  if (SD1==0) return SD2;
  while(cell_next(ptr)!=0) ptr=cell_next(ptr);
  cell_next(ptr)=SD2;
  return SD1;
}

/* end cly_cells.c */

/************************************************************************/
/***************** implementation code from cly_initial.c ***************/
/************************************************************************/

/*******************************************************************
** Initial Simplex Program
*******************************************************************/
/*
**  cell  Initial_simplex(node *S1, node *S2){
**
**  input:
**       A list of points S1.
**
**  output:
**       NULL if "PC" lies in a hyperplane.
**       A full dimensional simplex (S) from "PC" otherwise.
**
**  effects:
**         "B" points to a list containing all points of "PC"
**             wich are not in (S).
**
** method:
**   a) initialize (S) with  upper and lower hulls of PC with
**      respect to the first standard basis vector ei.
**   b) repeatedly  calculate a normal (N) to the span of (S),
**      and add one of the  extreme points of PC (with respect
**       to (N)) which does not already lie in the span of (S).
*/
static int check_normal(Ipnt,HMatrix,Ipnt *,Ipnt *);
#define CU cell_U(initial)
#define CM cell_H(initial)
#define CL cell_norm(initial)

cell  cly_initial_splx(Ipnt *S1, Ipnt *S2){
  cell initial;
  int i,r,j;
  Ipnt ptr=*S1,refpt=0,newpt,tmp;

  initial=cell_new(cly_N,cly_R);

  /*
  ** Initialize global matrices
  **   M will represent the linear space spanned by points
  **   allready chosen. (this will be an invarient of the loop).
  **   L will hold the next search direction, originally (0,...,0,1)
  */
  HMresize(CU,1,1);
  HMresize(CL,1,cly_N);
  HMresize(CM,1,cly_N);
  HLVset(CL,1,1);
  for(i=2;i<=cly_N;i++) HLVset(CL,i,0);
  /*
  ** find two points not normal to CL=(0,....,0,1); and update M
  */
  if (check_normal(ptr,CL,&refpt,&newpt)==FALSE) return FALSE;
  cell_pnt(initial,0)=refpt;
  cell_pnt(initial,1)=newpt;
  for(i=1;i<=cly_N;i++)
      HLMset(CM,1,i,cell_point(initial,1,i)-cell_point(initial,0,i));
  r=1;
  while (++r<=cly_N){
    /*
    ** Calculate new search direction L
    ** (normal to affine span of points allready chosen)
    */
    
    HMfactor(CM,CU);
    HMbacksolve(CM,CL);
    HMgcd_reduce(CL); 


    /*
    ** if possible pick a new point not in the same hyperplane,
    ** normal to L, as points allready chosen.
    */
    if (check_normal(ptr,CL,&refpt,&newpt)==FALSE){
 #ifdef LOG_PRINT     
      fprintf(stderr /* was cly_err */,"Failure in initial_simplex():");
      fprintf(stderr /* was cly_err */,"Point Config is not full dimensional\n")
#endif
;
       return 0;
    }
    cell_pnt(initial,r)=newpt;
    CM=HMsubmat(CM,r,cly_N);
    CU=HMresize(CU,r,r);
    for(i=1;i<=cly_N;i++)
       HLMset(CM,r,i,cell_point(initial,r,i)-cell_point(initial,0,i));
  }
  /* 
  ** reset sizes for matrices in cell 
  */
  HMresize(cell_norm(initial),1,cly_N+1);
  HMresize(cell_U(initial),cly_N,cly_N);
  HMresize(cell_H(initial),cly_N,cly_N+1);
  /*
  ** set normal to initial simplex
  **   (it is unlifted so norm = (0,...,0,1)
  */
  for(i=1;i<=cly_N;i++) HLVset(cell_norm(initial),i,0);
  HLVset(cell_norm(initial),cly_N+1,1);

  /* set up type vector */
  for(i=1;i<=cly_R;i++) cell_type(initial,i)=-1;
  for(i=0;i<=cly_N;i++) cell_type(initial,Ipnt_idx(cell_pnt(initial,i)))++;
  /*
  ** setup factorization matrix CU[v1-v0,....,vn-v0,0]=H
  */
  for(i=1;i<=cly_N;i++){
     for(j=1;j<=cly_N;j++){
        HLMset(cell_H(initial),j,i,
              cell_point(initial,i,j)-cell_point(initial,0,j));
     }
     HLMset(cell_H(initial),i,cly_N+1,0);
  }