pelclyal.cc

Gambit 是一个游戏库理论软件

  HMfactor(cell_H(initial),cell_U(initial));
  if (HLeq(HMget(cell_H(initial),cly_N,cly_N),0)){
            bad_error("initial simplex not full d");
  }
  cell_set_volume(initial);
  /*
  ** Set S2=points in initial simplex.
  **     S1=S1/S2.
  */
  ptr=*S1;
  *S1=0;
  *S2=0;
  while(ptr!=0){
     int In=FALSE;
     for(i=0;i<=cly_N;i++){
        if(ptr==cell_pnt(initial,i)) In=TRUE;
     }
     tmp=ptr;
     ptr=Ipnt_next(ptr);
     if (In==FALSE) {
         Ipnt_next(tmp)=*S1;
         *S1=tmp;
     }
     else{
         Ipnt_next(tmp)=*S2;
         *S2=tmp;
     }
  }
 return initial;
}


/*
** int check_normal(node PC,Imatrix norm, node *refpt, node *newpt){
**
** (auxilary function for initial_simplex)                        
**                                                                
**  input:                                                        
**    A point configuration "PC"                                  
**    A normal direction "N" and a reference point "refpt"        
**          (if "refpt" is null, "refpt" gets initalized to point 
**            to a point on the lower hull)                       
**          ("norm" and "refpt" together determin a hyperplane H) 
**                                                                
**  output:                                                       
**    FALSE if the point configuration lies hyperplane normal to N
**    TRUE otherwise.                                             
**                                                                
**  side effects:                                                 
**   "newpt" holds a point which has maximal distance from (H).   
**    if refpt was null it gets initialized (see "input").        
**                                                                
**  method:  just computes dot products and save a point acheiving
**           maximum and minimum values.                       
**           then compair these against value for refpt.
*/

static int check_normal(Ipnt PC,
                        HMatrix norm,
                        Ipnt *refpt, Ipnt *newpt){
  Hint dot,dmin,dmax;
  int i,first=TRUE;
  Ipnt pmin = NULL;
  Ipnt pmax = NULL;
  Hinit(dot,0);
  Hinit(dmin,0);
  Hinit(dmax,0);
  while(PC!=0){
    Hnorm_dot(dot,norm,PC);
    if (first==TRUE) {
       HHset(dmin,dot);
       HHset(dmax,dot);
       pmin=pmax=PC;
       first=FALSE;
    }
    else {
      if (HHlt(dot,dmin)){
          HHset(dmin,dot);
          pmin=PC;
      }
      else if (HHlt(dmax,dot)){
        HHset(dmax,dot);
        pmax=PC;
      }
    }
  PC=Ipnt_next(PC);
  }
  if (HHeq(dmin,dmax)) return FALSE;
  if (*refpt==0) {
     *refpt=pmin;
     *newpt=pmax;
  }
  else{
      Hnorm_dot(dot,norm,*refpt);
      if (HHeq(dot,dmin)) *newpt=pmax;
       else *newpt=pmin;
  }
  Hfree(dot);
  Hfree(dmin);
  Hfree(dmax);
  return TRUE;
}

/* end cly_initial.c */

/************************************************************************/
/***************** implementation code from cly_update.c ****************/
/************************************************************************/

static int is_flipped(cell c1, cell c2, int *i1, int *i2); 
static cell cell_pivot(cell C, Ipnt x, int idx);
static cell subdiv_add_cell(cell c, cell SDx);
int cell_find_lift(cell c, Ipnt pt);


/*
** is_flipped
**   Input:  cells   c1, and c2.
**   Output: TRUE if c1 and c2 are related by pivoting
**           FALSE otherwise.
**  Side Effects: i1 and i2 contain the indices of the pivoting
**                vertices for c1 and c2 respectivly.
*/
static int is_flipped(cell c1, cell c2, int *i1, int *i2){
  int i, found1=FALSE, found2=FALSE;
  for(i=0;i<=cly_N;i++){
    if (point_is_in(cell_pnt(c1,i),c2)==FALSE){
      if (found1==FALSE){
         *i1=i;
         found1=TRUE;
      }
      else return FALSE;
    }
    if (point_is_in(cell_pnt(c2,i),c1)==FALSE){
      if (found2==FALSE){
         *i2=i;
         found2=TRUE;
      }
      else return FALSE;
    }
  }
  if (found1==FALSE||found2==FALSE)
     bad_error("is_flipped: cell differences not detected");
  return TRUE;
}

/*
**   cell_pivot
**     Input.  A cell C.
**             A point x.
**             An index idx.
**    Output.  A copy of C, with point x substituted for vertex idx.
**             (the pointer field is set to indicate that cell C
**              and the new cell are related by a flip).
** 
*/
static cell cell_pivot(cell C, Ipnt x, int idx){
  int i,j;
  cell ncell;
  if (C==0) bad_error("null cell in cell_pivot\n");

  /* copy old cell */
  ncell=cell_new(cly_N,cly_R);
  for(i=1;i<=cly_R;i++) cell_type(ncell,i)=cell_type(C,i);
  for(i=0;i<=cly_N;i++) cell_pnt(ncell,i)=cell_pnt(C,i);

  cell_type(ncell,Ipnt_idx(cell_pnt(C,idx)))--;
  cell_type(ncell,Ipnt_idx(x))++;
  cell_pnt(ncell,idx)=x; /* substitute point x for vertex idx*/
  cell_ptr(ncell,idx)=C;   /* pivoting ncell at idx gives c*/
  cell_ptr(C,idx)=ncell;   /* pivoting c at idx gives ncell*/
  cell_fptr(ncell)--;
  cell_fptr(C)--;

  /* set up to calculate normal */
  /* load coordinate matrix */
  for(i=1;i<=cly_N;i++){
    for(j=1;j<=cly_N;j++){
      HLMset(cell_H(ncell),i,j,
             cell_point(ncell,i,j)-cell_point(ncell,0,j));
      }
    HLMset(cell_H(ncell),i,cly_N+1,cell_lift(ncell,i)-cell_lift(ncell,0));
  }
  
  /* calculate normal and remove common factors from coordinates*/
  HMfactor(cell_H(ncell),cell_U(ncell));
  HMbacksolve(cell_H(ncell),cell_norm(ncell));
  HMgcd_reduce(cell_norm(ncell));

  /* switch direction if necessary to ensure inner normal */
  if (HLlt(HVget(cell_norm(ncell),cly_N+1),0)) {
     for(j=1;j<=cly_N+1;j++) Hneg(HVget(cell_norm(ncell),j));
  }
  else if (HLeq(HVget(cell_norm(ncell),cly_N+1),0)) 
      bad_error("Cell perpendicular in Cell_norm");

 /* load and factor point matrix */
  for(i=1;i<=cly_N;i++){
     for(j=1;j<=cly_N;j++){
        HLMset(cell_H(ncell),j,i,
               cell_point(ncell,i,j)-cell_point(ncell,0,j));
     }
     HLMset(cell_H(ncell),i,cly_N+1,0);
  }
  HMfactor(cell_H(ncell),cell_U(ncell));

  cell_set_volume(ncell);
  
  return ncell;
}

/*
**  subdiv_add_cell
**    Input:     A cell c
**               A list of cells SDx
**               (with the invariant condition that all cells
**                of SDx related by a flip point to each other
**                through the relevent pointer fields.)
**   Output:     c appended to SDx
**               (with the invariant condition maintained).
**
**
*/
static cell subdiv_add_cell(cell c, cell SDx){
  int i,j;
  cell ptr;
  if (c==0) bad_error("null cell in subdiv_add_cell\n");
  if (SDx==0) return c;
  ptr=SDx;
  while(ptr!=0){
      if (is_flipped(c,ptr,&i,&j)==TRUE){
          if (cell_ptr(c,i)!=0)
              bad_error("new cell already flipped in sudiv_add_cell");
          if (cell_ptr(ptr,j)!=0)
              bad_error("old cell already flipped in sudiv_add_cell");
          cell_ptr(c,i)=ptr;
          cell_fptr(c)--;
          cell_ptr(ptr,j)=c;
          cell_fptr(ptr)--;
      }
      if (cell_next(ptr)==0){
              cell_next(ptr)=c;
              cell_next(c)=0;
              ptr=0;
      }
      else ptr=cell_next(ptr);
  }
  return SDx;
}

/*
**    New_Facets
**      Input        a subdivision SD
**                   a point x
**      Output       The new cells resulting by "placeing"
**                   x (implicitly using a conservitive lifting)
**
**      Method:        (let n1=n+1).
**             Given a cell C and an affine relation
**                  x=g0*v_0+...gn1*v_n1,     g1+...+gn1=1;
**
**             Lemma 2.20 of Verschelde et. al shows that
**             should be added by pivoting x into C in position
**             i iff gi<0.
**
**          For each Cell C.
**            The matrix M=(v_1-v_0,...,v_{n-1}-v_0,x-v_0) is formed
**            and the system M*L=0 is solved,
**            (giving g_i=-L_i/L_n and g_0=(L_0+...+L_n1)/L_n1
**
**            if (g_i < 0 and C has not allready been pivoted at v_i)
**                         (i.e. facet oposite v_i is an outer facet)
**            then
**             C is pivoted at v_i and vi is added to teh list of new
**             cells
**
**
*/
 cell cly_new_cells(cell D,cell Dl, Ipnt x){
  cell SDx=0,B;
  int i,j,n,tog;
  Hint l0,tmp;

  n=cly_N;
  Hinit(l0,0);
  Hinit(tmp,0);
  HHset(tmp,l0);
  HMresize(cly_L,1,cly_N+1);
  if (D==0){
    tog=1;
    B=Dl;
  }
  else {
    tog=0;
    B=D;
  }
  while(B!=0){
    if (cell_is_outer(B)==TRUE){
      /* load matrix L=(x-c0) */
      for(j=1;j<=n;j++)
        HLMset(cly_L,1,j,Ipnt_coord(x,j)-cell_point(B,0,j));
      /* put U*L in last collumb of H */
      for(i=1;i<=n;i++){
        HLMset(cell_H(B),i,cly_N+1,0);
        for(j=1;j<=n;j++){ 
          HHmul(tmp,HMget(cell_U(B),i,j),HMget(cly_L,1,j));
          HHadd(HMget(cell_H(B),i,cly_N+1),HMget(cell_H(B),i,cly_N+1),tmp);
        }   
      }
      /* solve H*L=0 */
      HMbacksolve(cell_H(B),cly_L);
      if (HLeq(HVget(cly_L,n+1),0)){
         cell_print(B);
         bad_error("simplex not full dim in new_facets()");
      }
      /*
      ** we have L1*(c1-c0)+...Ln*(cn-c0)+Ln1*(x-c0)
      **   that is we have l1c0+...+lncn=x for
      **     l0=(L1+...+Ln+Ln1)/Ln1;
      **     li=-Li/Ln1
      **    we are supposed to flip if li<0 so for i in 1...n
      **        we test sign(Li)*sign(Ln1)>0
      **    and for i=0 we test sign(Ln1)*sign(L1+...+Ln)<0
      */

      HHset(l0,HVget(cly_L,n+1));
      for(i=1;i<=n;i++){
       HHmul(tmp,HVget(cly_L,cly_N+1),HVget(cly_L,i));
       if ((HLgt(tmp,0)) && (cell_ptr(B,i)==0)){
          SDx=subdiv_add_cell(cell_pivot(B,x,i),SDx);
        }
        HHadd(l0,l0,HVget(cly_L,i));
      }
      HHmul(tmp,HVget(cly_L,cly_N+1),l0);
      if ((HLlt(tmp,0)) && (cell_ptr(B,0)==0)){
        SDx=subdiv_add_cell(cell_pivot(B,x,0),SDx);
      }
    }
    B=cell_next(B);
    if (B==0 && tog ==0){
      B=Dl;
      tog=1;
    }
  }
  Hfree(tmp);
  Hfree(l0);
  return SDx;
}

/*
**  cell_find_lift
**   Input: a cell c
**          a point pt
**   Output: lowest lifting value wich can be given s while
**           keeping it in upper half space defined by lifted cell
**           c.
**
**   Note:  cell_find_lift requires normals to be stored with cell.
*/
int cell_find_lift(cell c, Ipnt pt){
  Hint dotp,dotc,res;
  int i;
  HMatrix norm;
  Hinit(dotp,0);
  Hinit(dotc,0);
  Hinit(res,0);
  norm=cell_norm(c);
  Hnorm_dot(dotp,norm,pt);
  Hnorm_dot(dotc,norm,cell_pnt(c,0));
  HLmul(res,HVget(norm,cly_N+1),cell_lift(c,0));
  HHadd(dotc,dotc,res);
  HHsub(res,dotc,dotp);
  HHdiv(dotp,res,HVget(norm,cly_N+1));
  HLadd(res,dotp,1);
  Hfree(dotp);
  Hfree(dotc);
  if (HLlt(res,1)) return 1;
  return HtoL(res);
}
/*
** subdiv_find_lift
**   Input:   a subdivision s
**            a point pt
**  Output:   maximum lifting values pt can be given while keeping
**            it simultaniously in upper (open) half spaces defined
**            by cells of s.
*/
 int cly_find_lift(cell s, Ipnt pt){
 int res, temp;
 if (s==0) return 1;
 res=cell_find_lift(s,pt);
 while((s=cell_next(s))!=0)
  if (res<(temp=cell_find_lift(s,pt))) res=temp;
 return res;
}

/* end cly_update.c */

/************************************************************************/
/************** implementation code from cly_triangulate.c **************/
/************************************************************************/

Ipnt Internalize_Aset(aset A);
void free_globals(void);
static cell new_cayley_triangulate(Ipnt *PC);
static node subdiv_get_norms(cell S, Imatrix T, int *mv);

node cly_triangulate(aset A, Imatrix T, int order, int lift){
   int mv;
   Ipnt Pts;
   cell D;
   node Res=0;
   LOCS(2);
   PUSH_LOC(A);
   PUSH_LOC(Res);
   Pts=Internalize_Aset(A);
   cly_order = order;
   cly_lift = lift; 

   D=new_cayley_triangulate(&Pts);

/*   Ipnt_list_print(Pts);*/
   if (order==TRUE) order_subdiv(&D);
   subdiv_print(D);
   if (order==TRUE) print_all_volumes(D);
   if (lift==TRUE){
     lift_original_points(Pts);
     if (T!=0) Res=subdiv_get_norms(D,T,&mv);
     else Res=0;
   }
   else {
    warning("finding facets not supported\n");
 /*   Res=subdiv_to_facets(D);*/
   }
   subdiv_free(D);
   points_free(Pts);
   free_globals();
   POP_LOCS();
   return Res;
}

Ipnt Internalize_Aset(aset A){
  Ipnt point=0,points=0;
  node ptr,ptc,ptp;
  int r;
  cly_Dim=aset_dim(A)-1;
  cly_R=aset_r(A);
  cly_N=cly_Dim+cly_R-1;
  cly_U=HMnew(cly_N,cly_N);
  cly_M=HMnew(cly_N+1,cly_N+1);
  cly_L=HVnew(cly_N+1);
  Hinit(cly_temp,0);

  /* internalize points of Aset */
  ptr = aset_start_cfg(A);
  r=0;
  cly_Npts=0;
  while ((ptc = aset_next_cfg(&ptr)) != 0) {
      r++;
      while ((ptp = aset_next_pnt(&ptc)) != 0) {
        point=Ipnt_new(ptp,r);
        Ipnt_next(point)=points;
        cly_Npts++;
        points=point;
      }
  }
  return points;
}