pelpsys.cc

Gambit 是一个游戏库理论软件

/* psys.c */

/*
  This file now houses all the implementation code from files beginning with 
psys in the original distribution.
*/

#include "pelpsys.h"

/**************************************************************************/
/****************** implementation code from psys_aset.c ******************/
/**************************************************************************/

aset psys_to_aset(psys P){
  aset res=0,pnt=0;
  int r=1,j,d,t;
  char  *s,c='a'-1;
  LOCS(2);

  PUSH_LOC(res);
  PUSH_LOC(pnt);
  d=psys_d(P);
  r=psys_r(P);

  res=aset_new(r,d+1);

  FORALL_BLOCK(P,
    psys_Bstart_poly(P,psys_bkno(P));
    t=psys_eq_size(P);
    c++;
    FORALL_MONO(P,
      s=(char *)mem_malloc(6*sizeof(char));
      s[0]=c;
      sprintf(s+1,"%d",t--);
      pnt=aset_new_pt(d+1,s);
      for(j=1;j<=d;j++) aset_pnt_set(pnt,j,*psys_exp(P,j));
      aset_pnt_set(pnt,d+1,*psys_def(P));
      aset_add(res,psys_bkno(P),pnt);
    )
  )  
     
  POP_LOCS();
  return res;
}    


 
#define T(i) (*IVref(T,i))
#define Pcoord(i) (*IVref(pnt_coords(P),i))
psys aset_to_psys(aset A,Ivector T, int seed){
 int i,j,n=0,r=0;
 int R,N,M=0,eqno=0,blkno=0;
 psys Sys;
 node Aptr,C,Cptr,P;
 double t;
 
N=aset_dim(A)-1;

 if ((R=aset_r(A))!=IVlength(T))
      bad_error("Gen_Poly: Aset and Type Vector incompatible");
 
  rand_seed(seed);
 
  /* determine dimensions of psys */
  Aptr=aset_start_cfg(A);
  while((C=aset_next_cfg(&Aptr))!=0){
          r++;
          n+=T(r);
          M+=T(r)*pcfg_npts(C);
  }
  if (n!=N) bad_error("Gen_Poly: bad Type Vector");
  Sys=psys_new(N,M,R);
 
  /* set up semi-mixed structure data */
  Aptr=aset_start_cfg(A);
  while((C=aset_next_cfg(&Aptr))!=0){
    blkno++;
    *psys_block_start(Sys,blkno)=eqno;
    for(i=1;i<=T(blkno);i++){
      Cptr=aset_start_pnt(C);
      P=Car(Cptr);
      eqno++;
      r=0;
      do{ 
        psys_init_mon(Sys);
        t=rand_double(0,1);
        *psys_coef_real(Sys)=cos(2*PI*t);
        *psys_coef_imag(Sys)=sin(2*PI*t);
        *psys_def(Sys)=Pcoord(N+1);
        n=0;
        for(j=1;j<=N;j++){
           *psys_exp(Sys,j)=Pcoord(j);
           n+=*psys_exp(Sys,j);
        }
        if (r<n) r=n;/* calculate actual degree of current poly */
        *psys_homog(Sys)=-n; /* save degree of current monomial */
        psys_save_mon(Sys,eqno);
      }while((P=aset_next_pnt(&Cptr))!=0);
      /* Homogenize */
      psys_set_eqno(Sys,eqno);
      FORALL_MONO(Sys, (*psys_homog(Sys))+=r;)
    }
  }  
  return Sys;
}
#undef Pcoord
#undef T

/* end psys_aset.c */

/**************************************************************************/
/**************** implementation code from psys_binsolve.c ****************/
/**************************************************************************/

/*
**    copyright (c) 1995  Birk Huber
*/

/*
** psys_binsolve.c                             Birkett Huber 
**                                     Created On: 4-2-95
**                                  Last Modified: 8-6-95         
** Use Hermite normal form algorithm and gaussian elimination to 
** solve a simplicial system -- used to find initial solutions in 
** lifting homotopy.
*/

#define IMATRIX_FAST 1

#define CV_Set(V,i,j,re,im) {DMref((V),i,2*(j)-1)=(re);\
                             DMref((V),i,2*(j)  )=(im);}
#define CV_Read(V,i,j)   Complex(DMref(V,i,2*(j)-1),DMref(V,i,2*(j)))
#define CM_Set(M,i,j,re,im) {\
                DMref((M),2*(i)-1,2*(j)-1)=(re);\
                DMref((M),2*(i)  ,2*(j)  )=(re);\
                DMref((M),2*(i)-1,2*(j)  )=-1.0*(im);\
                DMref((M),2*(i)  ,2*(j)-1)=(im);}
#define S(i,j)  (*IMref(S,i,j))
#define U(i,j)  (*IMref(U,i,j))
#define Idx(i)  (*IVref(Idx,i))

/*
** psys_binsolve 
**       Input: a  polynomial vector P.
**      Output: On success a matrix of solutions.
**              On failure the null pointer (and an error message)
**              
*/
node psys_binsolve(psys sys){
 int n,i,j,k;
 int bstart=0,block_size;
 int m=1,row=0,cid;
 xpnt Dsol;
 Dmatrix Coef, Cf, Sol=0;
 Imatrix S,Idx,U=0;
 fcomplex ctmp;
 node Sol_List=0;
 LOCS(1);
 PUSH_LOC(Sol_List);

 n=psys_d(sys);

 /* 
 ** reserve and initialize auxilary space 
 */     
 Coef=Dmatrix_new(2*n,2*n);   /* non-constant monomial coefs  */
 Cf  =Dvector_new(2*n);     /* constant term for each equation */
 S   =Imatrix_new(n,n);   /* non-constant monomials exps in rows */
 U   =Imatrix_new(n,n); /* unitary matrix putting S into HNF */
 Idx =Ivector_new(n); /*controlls enumeration of ith roots of unity */

/*
** initialize coeficient matrices to zero
*/
 for(i=1;i<=2*n;i++){
  for(j=1;j<=2*n;j++){
      DMref(Coef,i,j)=0.0;
  }
  DVref(Cf,i)=0.0;
 }

/*
** Load matrix representation of Equation Coef*X^S = C 
*/
 psys_start_block(sys);
 psys_start_poly(sys);
 do {
   block_size=psys_block_size(sys);
   for(i=1;i<=block_size;i++){
     if (psys_start_mon(sys)==FALSE){
         warning("sys not simplicial in binsys (1)");
         goto cleanup;
     }
     if (i==1) for(k=1;k<=n;k++) Idx(k)=*psys_exp(sys,k);
     CV_Set(Cf,1,bstart+i,-1.0*(*psys_coef_real(sys)),
                        -1.0*(*psys_coef_imag(sys)))
     for(j=1;j<=block_size;j++){
        if (psys_next_mon(sys)==FALSE){
             warning("sys not simplicial in binsys (2)");
             goto cleanup;
        }
        CM_Set(Coef,bstart+i,bstart+j,*psys_coef_real(sys),
                                      *psys_coef_imag(sys))
        if (i==1){
           for(k=1;k<=n;k++){
              S(bstart+j,k)=*psys_exp(sys,k)-Idx(k);
           }
        }
     }
     psys_next_poly(sys);
   }
   bstart+=block_size;
 }while (psys_next_block(sys)==TRUE);

   

 /*
 ** find equivalent binomial sytem X^S=Cf
 */
 Dmatrix_Solve(Coef,Cf,2*n);

 /* 
 ** use hermite normal form to get a triangulated system X^(US)=C^U
 **   (use first row of Coef as temporary storage for calculating C^U)
 */
 for(i=1;i<=n;i++){
   DMref(Coef,1,2*i-1)=DVref(Cf,2*i-1);
   DVref(Cf,2*i-1)=1.0;
   DMref(Coef,1,2*i  )=DVref(Cf,2*i  );
   DVref(Cf,2*i  )=0.0;
 }
 /* replace S by its hermite normal form U*S (upper triangular)*/
 Imatrix_hermite(S,U);
 /* calculate C^U*/
 for(i=1;i<=n;i++){
   for(j=1;j<=n;j++){
     ctmp=Cmul(CV_Read(Cf,1,i),Cpow(CV_Read(Coef,1,j),U(i,j)));
     CV_Set(Cf,1,i,ctmp.r,ctmp.i)
   }
 }

 /* 
 ** Calculate number of solutions to be found 
 ** (determinant of triangular matrix S)
 */
 for(i=1;i<=n;i++){
   m*=(S(i,i));
   Idx(i)=1;
 }

 /*
 ** Reserve space for solutions and initialize them to vector of constant terms
 */
 Sol=Dmatrix_new(m,2*n);
 for(i=1;i<=m;i++){
   for(j=1;j<=2*n;j++){
     DMref(Sol,i,j)=DVref(Cf,j);
   }
 }

 /*
 ** Iterate through all the appropriate roots of unity while reverse solving
 */
 while(++row<=m){
   for(j=n;j>=1;j--){
      ctmp=Cmul(Croot(CV_Read(Sol,row,j),S(j,j)),
                    RootOfOne(Idx(j),S(j,j)));
      CV_Set(Sol,row,j,ctmp.r,ctmp.i)
      for(i=j-1;i>=1;i--){
         ctmp=Cmul(CV_Read(Sol,row,i),
                Cpow(CV_Read(Sol,row,j),-1*(S(i,j))));
         CV_Set(Sol,row,i,ctmp.r,ctmp.i)
      }
   }
   cid=n;
   while(cid<=n && cid >=1){
      if (Idx(cid)<S(cid,cid)){
        (Idx(cid))++;
        cid++;
      }
      else {
        Idx(cid)=0;
         cid--;
      }
   }
   if(cid<1 && row <m) warning("binsys stoping too soon");
 }

 /*
 ** Convert complex n-vector solutions to real 2-nvectors with
 ** and two extra coordinates for the implicit(complex) homogenizing
 ** parameter. and one coordinate for the deformation parameter
 */

 for(i=1;i<=m;i++){
    Dsol=xpnt_new(psys_d(sys));
    /* homog parameter */
    xpnt_h_set(Dsol,Complex(1.0,0.0));
    /* affine coords */
    for(j=1;j<=n;j++){
      xpnt_xi_set(Dsol,j,CV_Read(Sol,i,j));
    }
   /*deformation parameter*/
   xpnt_t_set(Dsol,0.0);
   Sol_List=Cons(atom_new((char *)Dsol,DMTX),Sol_List);
 }

/* 
** clean up and leave
*/
cleanup:
  Imatrix_free(Idx);
  Imatrix_free(S);
  Imatrix_free(U);
  Dmatrix_free(Cf);
  Dmatrix_free(Coef);
  Dmatrix_free(Sol);
  
POP_LOCS();
return Sol_List;
}
#undef CV_Set
#undef CV_Read
#undef CM_Set
#undef S
#undef U
#undef Idx

/* end psys_binsolve.c */

/**************************************************************************/
/******************* implementation code from psys_cont.c ******************/
/**************************************************************************/

/*******************************************************************
**  
**  Projective Newton Path Following                    Birk Huber
**                                                      Dec 31 1995
**  This is based on a description by T.Y. Li given in a lecture
**  at the SEA-95  (Nov 23-17) at CIRM
**
**  This version is written to allow for quick change to run as 
**  C-code. (i.e. only access/change singel elts in arrays).
**  uses auxilary functions norm, and GQR which must be modified to
**  take a basis indicating which collumns to use.
**  
** Changes
**   1: Endgame strategy  use dt=(1-t) if t+dt is > final tol
**      sugested strategy seemed to prevent actual completion.
**   2:(planned) try using a history mechanism so that some number
**     of consecutive good steps are required before increasing
**     dt. (it seems we take a lot of bad steps)
**   3:(planned)use an estimate on the condition number to decide
**     when to abort a path which is becomming singular.
**   4:(possible) put some mechanism to stop infinite loops in
**     the Newton iteration.
**   4:(possible) put in a mechanism to detect and abort paths
**     going to infinity.
*****************************************************************/


/* NOW REDUNDANT
#include <math.h>
#include "psys.h"
*/

/*
**     Controll Parameters
**        
**       PN_dt0         initial stepsize
**       PN_maxdt       maximum stepsize
**       PN_mindt       minimum stepsize
**       PN_scaledt     scaling factor for step size
**       PN_cfac        initial contraction required for approx root
**       PN_NYtol       Stop Newton when |F(X)|<Ntol
**       PN_NDtol       Stop final Newton when |DX|<PN_FDtol
**       PN_Nratio      Stop Newton when |DX|/|F(X)| < PN_Nratio
**       PN_tfinal      Destination Value. (just short of 1)
**       PN_FYtol       Stop final Newton when |F(X)|<PN_FYtol 
**       PN_FDtol       Stop final Newton when |DX|<PN_FDtol
**       PN_Fratio      Stop final Newton when |DX|/|F(X)| < PN_Fratio
**       PN_maxsteps    Give up after PN_maxsteps iterations of main loop
**        
*/

double PN_dt0=.01;        
double PN_maxdt=0.1;        
double PN_mindt=1E-14;  
double PN_scaledt=2;        
double PN_cfac=10;      
double PN_NYtol=1E-8;     
double PN_NDtol=1E-12;     
double PN_Nratio=1E-4;   
double PN_tfinal=1-1E-10;
double PN_FYtol=1E-10;   
double PN_FDtol=1E-14;   
double PN_Fratio = 1E-8; 
int    PN_maxsteps=1000; 
#define Zero_Tol 1E-13
/*
**  Auxilary Macroes/Functions
**
**  Backsolve(R,QY,DX,i,j)
**     R- NxN upper triangular, QY N-vector DX n-vector, i,j integers
**     solve  R*DX=QY  (using i,j as loop counters)
*/

#define Basis(i) (*IVref(Basis,(i)))
#define Z(i)     (*IVref(    Z,(i)))
#define X(i)     (DVref(    X,(i)))
#define Y(i)     (DVref(Y,i))
#define QY(i)    (DVref(QY,i))
#define DX(i)    (DVref(DX,i))
#define XOld(i)  (DVref(XOld,i))
#define R(i,j)   (DMref(R,i,j))
#define QT(i,j)  (DMref(QT,i,j))
#define JC(i,j)  (DMref(JC,i,j))
#define A(i,j)   (DMref(A,i,j))
#define Q(i,j)   (DMref(Q,i,j))


#define Backsolve(R,QY,DX,i,j) \
    for(i=N; i>=1; i--){ \
    DX(i)=0.0;\
    for(j=N;j>i;j--){DX(i)=DX(i)+R(i,j)*DX(j);}\
    DX(i)=(QY(i)-DX(i))/(R(i,i));\
  }

   
double norm(Dvector X){
  int i;
  double abs=0.0,tmp;
  for(i=1;i<=DVlength(X);i++){
        tmp=X(i);
        abs+=(tmp*tmp);
  }
  return sqrt(abs);
}

double BGQR(Dmatrix A,Ivector Basis,int N,Dmatrix Q,Dmatrix R){
  int i,j,k;
  double s,s1,s2,t1,t2, max_e=-1,min_e=-1;
  
  for (i=1;i<=N;i++){
    for (j=1;j<=N;j++) {
       Q(i,j)=0.0;
       R(i,j)=A(i,Basis(j));
    }
    Q(i,i)=1.0;
  }
  for(i=1;i<=N-1;i++){
    for(k=i+1;k<=N;k++){
      if (R(k,i)>Zero_Tol || R(k,i)<-Zero_Tol){
         s2=R(k,i);      s1=R(i,i);
         s=sqrt(s1*s1+s2*s2);
         s1=s1/s;        
         s2=s2/s;
         for(j=1;j<=N;j++){
           t1= s1*R(i,j)+s2*R(k,j);
           t2=-s2*R(i,j)+s1*R(k,j);
           R(i,j)=t1;
           R(k,j)=t2;
           t1= s1*Q(i,j)+s2*Q(k,j);
           t2=-s2*Q(i,j)+s1*Q(k,j);
           Q(i,j)=t1;
           Q(k,j)=t2;
         }
      }
    }
    s1=fabs(R(i,i));
    if (s1>max_e||max_e<0) max_e=s1;
    else if(s1<min_e || min_e<0) min_e=s1;
  }
  return max_e/min_e;
}                        
/*
**       Variable Definition and Initialization
*/

int psys_cont(Dvector X, psys P){
Dvector Y=0, QY=0,DX=0,XOld=0;
Dmatrix JC=0,QT=0,R=0;
Ivector Basis=0,Z=0;
int T;         /* total number of real variables appearing*/
int N;     /* total number of variables for 1 affine patch*/
int tsteps=0;
int i,j;
int max_idx,tmp;
int Rval=0;
double max_c,new_c;
double Nold,Nnew = 0.0;
double arclen=0.0,arcstep;
double dt=PN_dt0;
double Ndx=0.0;
double cond_num;

N=2*psys_d(P);
T=N+3;
if (DVlength(X)!=T) bad_error("Bad Starting point");

Z    =Ivector_new(2);
Basis=Ivector_new(N);
XOld =Dvector_new(T);
DX   =Dvector_new(N);
Y    =Dvector_new(N);
QY   =Dvector_new(N);
JC   =Dmatrix_new(N,T);
QT   =Dmatrix_new(N,N);
R    =Dmatrix_new(N,N);
Z(1)=1; Z(2)=2;
for(i=3;i<=T-1;i++) Basis(i-2)=i;

/*
**      Main Predictor Corrector Loop
*/
while (tsteps<PN_maxsteps && X(T)<=PN_tfinal){
  tsteps=tsteps+1;   

  /* 
  ** End game Step size (THIS DIFFERS FROM LI'S SUGGESTION)
  */
  if (X(T)+dt > PN_tfinal) dt=1-X(T);
  
  /*
  ** Choose Basis:
  ** A) Find complex coordinate with largest norm 
  */
  max_c = X(Z(1))*X(Z(1))+X(Z(2))*X(Z(2));
  max_idx=0;
  for(i=1;i<=N;i=i+2){
    new_c=X(Basis(i))*X(Basis(i))+X(Basis(i+1))*X(Basis(i+1));
    if (new_c>max_c){
      max_c=new_c;
      max_idx=i;
    }
  }
  /*
  ** B) If largest coord is not allready homog param make it so.
  */
  if (max_idx>0){