pelclyal.cc

Gambit 是一个游戏库理论软件

void free_globals(void){
/* should also reset all global params*/
 HMfree(cly_L);
 HMfree(cly_M);
 HMfree(cly_U);
 Hfree(cly_temp);
 next_id=1;
}

/*
**  Alg 2.10 of Vershelde
**
**  Input: A point configuration (with an extra coordinate to hold
**         lifting values (weights) should be initialized to zeros
**
**  Output: A triangulation of the subdivision (represented by
**          A list of point configurations)
*/
static cell new_cayley_triangulate(Ipnt *PC)
{
    Ipnt x=0,S1=0,S2=0;
    cell Dl=0;
 
    S1=*PC;
    if ((Dl=cly_initial_splx(&S1,&S2))==0)
         bad_error("failure in initial simplex");

    while(S1!=0){
      x=S1;
      S1=Ipnt_next(x);
      if (cly_lift==TRUE) Ipnt_lift(x)=cly_find_lift(Dl,x);          
      Dl=cly_subdiv_union(Dl,cly_new_cells(0,Dl,x));
      Ipnt_next(x)=S2;
      S2=x;
   }
   *PC=S2;
   return Dl;
}

/*
** subdiv_get_norms(cell S, Imatrix T, *int mv);
**
**     Input: a subdivision S (not sorted), and a type matrix T.
**    Output: a list of lifted inner normals to cells of S with
**            type T.
**    side effects: All space reserved for S is freed.
**                  mv holds the type T mixed volume.
*/
#define Norm(i) (*IVref(Norm,i))
static node subdiv_get_norms(cell S, Imatrix T, int *mv){
   node Nlist=0;
   Imatrix Norm;
   int i;
   LOCS(1);
   PUSH_LOC(Nlist);
   *mv=0;
   while (S!=0){
     if (Imatrix_order(T,cell_T(S))==0){
         Norm=Ivector_new(cly_Dim+1);
         for(i=1;i<=cly_Dim;i++) Norm(i)=HLVget(cell_norm(S),i);
         Norm(cly_Dim+1)=HLVget(cell_norm(S),cly_N+1);
         Nlist=Cons(atom_new((char *)Norm,IMTX),Nlist);
         *mv+=cell_volume(S);
     }
     S=cell_next(S);
   }
   POP_LOCS();
   return Nlist;
}

/* end cly_triangulate.c */

/************************************************************************/
/**************** implementation code from cly_continue.c ***************/
/************************************************************************/

static cell cayley_continue(Ipnt *PC, int threshold, int tweak);
static Ipnt Internalize_Aset_Cont(aset,Imatrix,int);
       Ipnt Internalize_Aset();
static void free_continuation_globals();
       void free_globals();
static void Update_Poly(Ipnt x);
static void Update_Solutions(cell Dl, int tweak);
static node relift(node *, int tweak);
static void norm_reset(cell ncell);
int cell_find_lift(cell,Ipnt);

static psys     Poly_Sys=0;
static node     Poly_Sols=0;
static Imatrix  Poly_Type=0;
static Imatrix  Poly_Norm=0;
static Imatrix  Poly_TNorm=0;
static Ipnt     *Poly_Pnts=0;
#define PPnts(i) (Poly_Pnts[(i)-1])
static HMatrix   Poly_H=0;
static HMatrix   Poly_U=0;
#define Poly_Type(i) (*IVref(Poly_Type,i))
#define Poly_Norm(i) (*IVref(Poly_Norm,i))
#define Poly_TNorm(i) (*IVref(Poly_TNorm,i))
/******************************************************************
**  Cayley Triangulation and Continuation
******************************************************************/

psys Cayley_continue(aset A,Imatrix T,node *Sols,int seed,int tweak){
   Ipnt S;
   cell D;
   psys res;
   S=Internalize_Aset_Cont(A,T,seed);
   D=cayley_continue(&S,10,tweak);
   subdiv_free(D);
   points_free(S);
   *Sols=Poly_Sols;
   res=Poly_Sys;
   free_continuation_globals();
   return res; 
}

#define T(i) (*IVref(T,i))
Ipnt Internalize_Aset_Cont(aset A, Imatrix T,int seed){
  Ipnt Pts, res;
  int i,tmp,max_monomials=0;
  /* 
  ** seed random number generator
  */
  rand_seed(seed);
  /* 
  ** initialize regular triangulation globals
  */
  res=(Pts=Internalize_Aset(A));
  Poly_Pnts=(Ipnt *)mem_malloc(cly_Npts*sizeof(Ipnt));
  /* 
  ** Store type vector in globally accessable location
  */
  Poly_Type=T;
  Poly_Norm=Ivector_new(cly_Dim+1);
  Poly_TNorm=Ivector_new(cly_Dim+1);
  Poly_H=HMnew(cly_N,cly_N+1);
  Poly_U=HMnew(cly_N,cly_N);
  /*
  ** Count total number of monomials needed
  */
  i=1;
  while(Pts!=0){
    PPnts(i++)=Pts;
    max_monomials+= T(Ipnt_idx(Pts));
    Pts=Ipnt_next(Pts);
  }
  /*
  ** initialize new psys- create it set up the block structure
  **                      implied by the type vector T.
  */
  Poly_Sys=psys_new(cly_Dim,max_monomials,cly_R);
  tmp=1;
  for(i=1;i<=IVlength(Poly_Type)+1;i++){
     *psys_block_start(Poly_Sys,i)=tmp;
     tmp+=Poly_Type(i);
  }
return res;
}
#undef T

void free_continuation_globals(){
  Poly_Sys=0;
  Poly_Type=0;
  Imatrix_free(Poly_Norm);
  Imatrix_free(Poly_TNorm);
  HMfree(Poly_U);
  HMfree(Poly_H);
  free_globals();
}


 

/*
**  Alg 2.10 of Vershelde
**
**  Input: A point configuration (with an extra coordinate to hold
**         lifting values (weights) should be initialized to zeros
**
**  Output: A triangulation of the subdivision (represented by
**          A list of point configurations)
*/
static cell cayley_continue(Ipnt *PC, int threshold, int tweak)
{
    Ipnt pt, x=0,S1=0,S2=0;
    cell D=0,Dl=0;
    int l;
    LOCS(1);
    PUSH_LOC(Poly_Sols);

    S1=*PC;
    D=cly_initial_splx(&S1,&S2);
    for(l=0;l<=cell_n(D);l++) Update_Poly(cell_pnt(D,l));
    
    while(S1!=0){
      x=S1;
      S1=Ipnt_next(x);
      Ipnt_lift(x)=(l=cly_find_lift(Dl,x));
      Update_Poly(x);
      if (l>threshold){
        Update_Solutions(Dl,tweak); 
 #ifdef LOG_PRINT
	fprintf(stdout /* was cly_out */,"flattening\n")
#endif
;
        for(pt=S2;pt!=0;pt=Ipnt_next(pt)) Ipnt_lift(pt)=0;
        psys_lift(Poly_Sys,0);
        D=cly_subdiv_union(D,Dl);
        Dl=0;
        Ipnt_lift(x)=1;
      }
      Dl=cly_subdiv_union(Dl,cly_new_cells(D,Dl,x));
      Ipnt_next(x)=S2;
      S2=x;
   }
   Update_Solutions(Dl,tweak);
   D=cly_subdiv_union(D,Dl);
   *PC=S2;
   POP_LOCS();
   return D;
}

  
/*
** Update_Poly(Ipnt x)  a monomail representing x is chosen with
**                      a random coefficient and added to the poly.
*/
static void Update_Poly(Ipnt x){
     int j,n;
     double t;
     
     n=psys_d(Poly_Sys);
     psys_Bstart_poly(Poly_Sys,Ipnt_idx(x));
     do{
       psys_init_mon(Poly_Sys);
       *psys_aux(Poly_Sys)=x;
       t=rand_double((int)0.0,(int)(2*PI));
       *psys_coef_real(Poly_Sys)=cos(t);
       *psys_coef_imag(Poly_Sys)=sin(t);
       *psys_def(Poly_Sys)=Ipnt_lift(x);
       for(j=1;j<=n;j++){
           *psys_exp(Poly_Sys,j)=Ipnt_coord(x,j);
       }
       psys_save_mon(Poly_Sys,psys_eqno(Poly_Sys));
     } 
     while(psys_Bnext_poly(Poly_Sys,Ipnt_idx(x))==TRUE);
}

int comp(node g1,node g2){
  return Imatrix_order(cell_norm((cell)Car(g1)),
                      cell_norm((cell)Car(g2)));
}

void update_list_insert(node g, node * L)
{
  node ptr = *L;
  LOCS(2);
  PUSH_LOC(*L);
  PUSH_LOC(g);
  if (g==0) bad_error("null node passed to update list");
  if (*L == 0 || comp(g, Car(*L)) > 0) *L = Cons(g, *L);
  else {
    while ((Cdr(ptr) != 0) && comp(g, Car(Cdr(ptr))) <= 0) 
       ptr = (node) Cdr(ptr);
    node_set_ptr(ptr,Cons(g, Cdr(ptr)), NODE, RIGHT);
  }
  POP_LOCS();
}

/*
** Udate_Solutions(cell D)
**                      use continuation on all solutions already 
**                      found.
**                      then for each cell in D do a lifting homotopy
**                      to add more solutions.
**                at end flatten all monomials, and solutions.
*/
static void Update_Solutions(cell Dl, int tweak){
   node tmp=0,lhead=0;
   LOCS(1);
   PUSH_LOC(lhead);
   tmp=Poly_Sols;
   /* update old solutions */
   while(tmp!=0){
       xpnt_t_set((Dmatrix)Car(Car(tmp)),0.0);
       tmp=Cdr(tmp);
   }

   if (Poly_Sols!=0){
  #ifdef LOG_PRINT     
     fprintf(stdout /* was cly_out */,"updating old cells\n")
#endif
;
       Poly_Sols=psys_hom(Poly_Sys,Poly_Sols,tweak); 
   }     

   /* find new solutions */

   /* order list of good cells by normal */
   /*      (with low numbered cells first) */
   while(Dl!=0){
     if (Imatrix_equal(Poly_Type,cell_T(Dl))==TRUE){
       update_list_insert(atom_new((char *)Dl,CELL),&lhead); 
     }
     Dl=cell_next(Dl);
   }
   /* cut out blocks, of cells */
   while(lhead!=0){
      Poly_Sols=list_cat(relift(&lhead,tweak),Poly_Sols);
   }   
   psys_lift(Poly_Sys,0);
   POP_LOCS();
}
   
/*
** 
*/
#define Tmp_Norm(i) (*IVref(Tmp_Norm,i))
static node relift(node *lhead, int tweak){
    Imatrix Tmp_Norm;
    int i,ht,tmp_ht;
    cell Tmp_Cell;
    psys Lead_Sys=0, Norm_Sys=0;
    node Lead_Sols=0,list2=0,lptr=0;
    Ipnt Mpt;
    LOCS(3);
    PUSH_LOC(*lhead);
    PUSH_LOC(Lead_Sols);
    PUSH_LOC(list2);
    /*
    ** check if cell list is empty
    */
    if (*lhead==0) bad_error("null list in relift");
    /*
    ** make copy of truncated original norm
    */
    Tmp_Norm=cell_norm((cell)Car(Car(*lhead)));
    for(i=1;i<=cly_Dim;i++){
       Poly_Norm(i)=Tmp_Norm(i);
    }
    Poly_Norm(cly_Dim+1)=Tmp_Norm(cly_N+1);
    /* 
    ** set up transformed and leading equations 
    */
    Norm_Sys=psys_norm_sub(psys_copy(Poly_Sys),Poly_Norm);
    Lead_Sys=psys_lead(Norm_Sys);
    /*
    ** initialize all points to unseen.
    */ 
    for(i=1;i<=cly_Npts;i++){
        Ipnt_seen(PPnts(i))=FALSE;
    }
    /*
    ** relift points, and cells
    */
    Tmp_Cell=(cell)Car(Car(*lhead));
    for(i=0;i<=cly_N;i++){
     Ipnt_seen(cell_pnt(Tmp_Cell,i))=TRUE;
     Ipnt_lift(cell_pnt(Tmp_Cell,i))=0;
    }
    while (Cdr(*lhead)!=0 && comp(Car(*lhead),Car(Cdr(*lhead)))==0){
      norm_reset(Tmp_Cell);
      list_push(list_pop(lhead),&list2);
      Tmp_Cell=(cell)Car(Car(*lhead));
      for(i=0;i<=cly_N;i++){
        if (Ipnt_seen(cell_pnt(Tmp_Cell,i))!=TRUE){
          lptr=list2;
          ht=cell_find_lift((cell)Car(Car(lptr))
                              ,cell_pnt(Tmp_Cell,i));
          while ((lptr=Cdr(lptr))!=0){
              tmp_ht=cell_find_lift((cell)Car(Car(lptr))
                                    ,cell_pnt(Tmp_Cell,i));
              if (ht<tmp_ht) ht=tmp_ht;
          }
          Ipnt_seen(cell_pnt(Tmp_Cell,i))=TRUE;
          Ipnt_lift(cell_pnt(Tmp_Cell,i))=ht;
        }
      }
    }
    norm_reset(Tmp_Cell);
    list_push(list_pop(lhead),&list2);

   /* lift monomials */
    FORALL_POLY(Lead_Sys,
      FORALL_MONO(Lead_Sys,
        Mpt=(Ipnt)(*psys_aux(Lead_Sys));
      if (Ipnt_seen(Mpt)==TRUE){
       *psys_def(Lead_Sys)=Ipnt_lift(Mpt);
      }
      else{
        lptr=list2;
        ht=cell_find_lift((cell)Car(Car(lptr)),Mpt);
        while ((lptr=Cdr(lptr))!=0){
           tmp_ht=cell_find_lift((cell)Car(Car(lptr)),Mpt);
           if (ht<tmp_ht) ht=tmp_ht;
        }
        Ipnt_seen(Mpt)=TRUE;
        Ipnt_lift(Mpt)=ht;
        *psys_def(Lead_Sys)=Ipnt_lift(Mpt);
      }
      );
    );
    /* 
    ** take solutions to Lead_Sys, as starting points for
    ** Norm_Sys and use homotopy.
    */
    lptr=list2;
    while (lptr!=0){
     Tmp_Norm=cell_norm((cell)Car(Car(lptr)));
     for(i=1;i<=cly_Dim;i++){
       Poly_TNorm(i)=Tmp_Norm(i);
     }
     Poly_TNorm(cly_Dim+1)=Tmp_Norm(cly_N+1);
#ifdef LOG_PRINT
     fprintf(stdout /* was cly_out */,"deforming to face:\n")
#endif
;
     Lead_Sols=list_cat(psys_solve(Lead_Sys,Poly_TNorm,tweak),Lead_Sols);
    lptr=Cdr(lptr);
    }
   /* 
   ** reset t value for these and use continuation on Sys_Norm,
   ** and return 
   */
   lptr=Lead_Sols;
   while(lptr!=0){
       xpnt_t_set((Dmatrix)Car(Car(lptr)),0.0);
       lptr=Cdr(lptr);
   }
 #ifdef LOG_PRINT
   fprintf(stdout /* was cly_out */,"deforming from face:\n")
#endif
;
   if (Lead_Sols!=0) Lead_Sols=psys_hom(Norm_Sys,Lead_Sols,tweak);
   psys_free(Lead_Sys);
   psys_free(Norm_Sys);
   POP_LOCS();
   return Lead_Sols;}
#undef Tmp_Norm

static void  norm_reset(cell ncell){
   int i,j;
  /* load coordinate matrix */
  for(i=1;i<=cly_N;i++){
    for(j=1;j<=cly_N;j++){
      HLMset(Poly_H,i,j,
             cell_point(ncell,i,j)-cell_point(ncell,0,j));
    }
 HLMset(Poly_H,i,cly_N+1,cell_lift(ncell,i)-cell_lift(ncell,0));
  }

  /* calculate normal and remove common factors from coordinates*/
  HMfactor(Poly_H,Poly_U);
  HMbacksolve(Poly_H,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  reset norm");
}
#undef Poly_Type
#undef Poly_Norm
#undef Poly_TNorm

/* end cly_continue.c */