pelpsys.cc

Gambit 是一个游戏库理论软件

    tmp=Basis(max_idx);
    Basis(max_idx)=Z(1);
    Z(1)=tmp;
    tmp=Basis(max_idx+1);
    Basis(max_idx+1)=Z(2);
    Z(2)=tmp;
  }

  /*
  **  Store Current Solution: (To allow reset if step is rejected)
  */
  for (i=1;i<=T;i++) XOld(i)=X(i);

  /*
  ** Predictor:  Euler predictor step 
  ** A) Solve system JC(:,Basis)DX=-JC(:,T)
  */
  JC=psys_jac(P,X,JC);
  cond_num=BGQR(JC,Basis,N,QT,R);
  /* 
  ** Form RHS 
  */
  for(i=1;i<=N;i++){
    QY(i)=0.0;
    for(j=1;j<=N;j++){
      QY(i)=QY(i)-QT(i,j)*JC(j,T);
    }
  }
  /* 
  ** Backsolve
  */
  Backsolve(R,QY,DX,i,j)

  /*
  ** B) incriment X(T), and X(Basis).
  */
  X(T)=X(T)+dt;
  for(i=1;i<=N;i++) X(Basis(i))=X(Basis(i))+dt*DX(i);

  /* Corrector First Step */
  Y=psys_eval(P,X,Y);
  Nold=norm(Y);             
  JC=psys_jac(P,X,JC); 
  cond_num=BGQR(JC,Basis,N,QT,R);
  /* form RHS: */
  for(i=1;i<=N;i++){
    QY(i)=0;
    for(j=1;j<=N;j++){
        QY(i)=QY(i)-QT(i,j)*Y(j);
    }
  }
  /* Back solve */
  Backsolve(R,QY,DX,i,j)

  /* update */
  for(i=1; i<=N;i++) X(Basis(i))=X(Basis(i))+DX(i);
  Y=psys_eval(P,X,Y); 
  Nnew=norm(Y);

  /*Decide if performance is o.k.*/
  if (Nnew>(Nold/PN_cfac)&&(Nnew>PN_NYtol)){
    /*if performance is not o.k. decrease step size and try again*/
    if (dt>PN_mindt){
      dt=dt/PN_scaledt;     
      for(i=1;i<=T;i++)X(i)=XOld(i);
    }
    else {
      fprintf(stdout /* was psys_logfile */,"Minimum Step size reached\n");
      printf("Minimum Step size reached\n");
      Rval=1;
      goto cleanup;
    }
    fprintf(stdout /* was psys_logfile */,"Stp=%d, T=%g, dt=%g\n", tsteps,X(T),dt);
  }
  else{     
    /* 
    ** if performance is o.k. run Corrector loop
    ** (Maybe a good idea to limit number of Newton Steps)
    */
    Ndx=norm(DX);
    while (Nnew>PN_NYtol && Ndx>PN_NDtol && Ndx/Nnew > PN_Nratio){
    JC=psys_jac(P,X,JC); 
    cond_num=BGQR(JC,Basis,N,QT,R);
      /* form RHS*/
      for(i=1;i<=N;i++){
        QY(i)=0;
        for(j=1;j<=N;j++) QY(i)=QY(i)-QT(i,j)*Y(j);
      }
      Backsolve(R,QY,DX,i,j)

      for (i=1;i<=N;i++) X(Basis(i))=X(Basis(i))+DX(i);
      Y=psys_eval(P,X,Y); 
      Nnew=norm(Y);
      Ndx=norm(DX);
    }

    /*add size of current step to arclenth approx*/
    arcstep=0.0;
    for(i=1;i<=T;i++) arcstep+=pow(X(i)-XOld(i),(long)2);
    arclen=arclen+sqrt(arcstep);

    /* increase step size if allowed*/
    if (dt <PN_maxdt) dt=dt*PN_scaledt;      
    fprintf(stdout /* was psys_logfile */,"Stp=%d, T=%g, Y=%g, arclen=%f\n", 
           tsteps,X(T),Nnew,arclen);
  }
}

if (tsteps>PN_maxsteps){
 printf(" too many steps taken\n");
 Rval=1;
 goto cleanup;
}
  
/* RUN NEWTON's METHOD TO FINAL TOLERANCE (IF DESIRED).*/
  X(T)=1.0;
  Y=psys_eval(P,X,Y); 
  Nnew=norm(Y);
  if (Nnew>PN_FYtol){
    fprintf(stdout /* was psys_logfile */,
           "Starting final Newton iteration: T=%g,Y=%g\n",X(T),Nnew);
   do {
     JC=psys_jac(P,X,JC); 
     cond_num=BGQR(JC,Basis,N,QT,R);
     /* form RHS: */
     for(i=1;i<=N;i++){
       QY(i)=0;
       for(j=1;j<=N;j++) QY(i)=QY(i)-QT(i,j)*Y(j);
      }
      Backsolve(R,QY,DX,i,j)
      for(i=1;i<=N;i++) X(Basis(i))=X(Basis(i))+DX(i);
      Y=psys_eval(P,X,Y); 
      Nnew=norm(Y);
      Ndx=norm(DX);
      fprintf(stdout /* was psys_logfile */,
        "     Y=%g, Ndx=%g, Ndx/Y=%g cond=%g\n",Nnew,Ndx,Ndx/Nnew,cond_num);
   }while (Nnew>PN_FYtol && Ndx >PN_FDtol && Ndx/Nnew > PN_Fratio);
  }

cleanup:
printf("Stp=%d, T=%g,Y=%g, arclen=%f, FLAG=%d\n", 
           tsteps,X(T),Nnew,arclen,Rval);
Ivector_free(Z);
Ivector_free(Basis);
Dvector_free(XOld);
Dvector_free(DX);
Dvector_free(Y); 
Dvector_free(QY);
Dmatrix_free(JC);
Dmatrix_free(QT);
Dmatrix_free(R);

return Rval;}


/* end psys_cont.c */

/**************************************************************************/
/******************* implementation code from psys_def.c ******************/
/**************************************************************************/

/*
** Psys class definition. 
**    -  square system of equations, listed by blocks
*/
struct psys_t{
    int N;     /* number of variables */
    int Neq;   /* number of equations -- usually Neq==N*/
    int R;     /* number of different support sets */
    int Mmax;  /* Max number of Monomials allowed */
    int M;     /* total number of monomials used */
    int *estart;
    int *bstart;
    int *degrees;
    int *tops;  /* a vector of info desribing system */ 
    int *istore; /* an Mx(N+3) matrix of exponents for monomials*/
 double *dstore; /* and Mx2 matrix of coeficients for monomials */
    int curr_eqn;
    int curr_mon;
    int curr_block;
    int free_mon;
    void **aux;
   Dvector Trans;
};
/* integer and double  and next pointer part of ith monomial*/
#define MonI(P,i) ((P)->istore+(((i)-1)*(((P)->N)+3)))
#define MonD(P,i) ((P)->dstore+((i)-1)*2)
/* access macroes for ith monomial */
#define Mon_coefR(P,i)   ((MonD(P,i))[0])
#define Mon_coefI(P,i)   ((MonD(P,i))[1])
#define Mon_homog(P,i)   ((MonI(P,i))[0])
#define Mon_defv(P,i)    ((MonI(P,i))[1])
#define Mon_next(P,i)    ((MonI(P,i))[2])
#define Mon_exp(P,i,d)   ((MonI(P,i))[2+(d)])
#define Mon_aux(P,i)     (((P)->aux)[i-1])
#define Mon_curr(P)      ((P)->curr_mon)

/* access macroes for ith equation */
#define Eqn_start(P,i)    (((P)->estart)[(i)-1])
#define Eqn_size(P,i) (((P)->tops)[(i)-1])
#define Eqn_curr(P)   ((P)->curr_eqn)
#define Eqn_deg(P,i)    ((P)->degrees[(i)-1])
/* access macroes for the ith block */
#define Blk_start(P,i)    (((P)->bstart)[(i)-1])
#define Blk_size(P,i)     (Blk_start(sys,i+1)-Blk_start(sys,i))
#define Blk_curr(P)       ((P)->curr_block)
/* access macroes for system params */
#define Sys_M(P)       ((P)->M)
#define Sys_Mmax(P)    ((P)->Mmax)
#define Sys_Neq(P)        ((P)->Neq)
#define Sys_N(P)          ((P)->N)
#define Sys_R(P)          ((P)->R)
#define Sys_Mon_New(P)    ((P)->free_mon)

/*
** constructor/destructor functions:
*/
#define MON_SIZE (m*(n+3+4))
#define ISTORE_SIZE  ((MON_SIZE+3*n+r+1))
psys psys_new(int n, int m, int r){
    int i;
    psys res;
    res=(psys)mem_malloc(sizeof(struct psys_t));
    res->istore=(int *)mem_malloc(ISTORE_SIZE*sizeof(int));
    res->dstore=(double *)mem_malloc(2*m*sizeof(double));
    res->aux=(void **)mem_malloc(m*sizeof(void *)); 
    Sys_N(res)=n;
    Sys_Neq(res)=n;
    Sys_R(res)=r;
    Sys_M(res)=0;
    Sys_Mmax(res)=m;
    res->estart=res->istore+MON_SIZE;
    res->tops=res->estart+n;
    res->degrees=res->tops+n;
    res->bstart=res->degrees+n;
    for(i=0;i<ISTORE_SIZE;i++)res->istore[i]=0;
    Blk_start(res,r+1)=n+1;
    /* put all monomials on free list */
    Sys_Mon_New(res)=1;
    for(i=1;i<=m;i++){
        Mon_next(res,i)=i+1;
        Mon_aux(res,i)=0; 
        Mon_coefR(res,i)=0.0;
        Mon_coefI(res,i)=0.0;
    }
    res->curr_mon=0;
    res->curr_eqn=0;
    res->curr_block=0;
    return res;
}

void psys_free(psys sys){
  mem_free((void *)(sys->dstore));
  mem_free((void *)(sys->istore));
  mem_free((void *)(sys->aux)); 
  mem_free((void *)(sys));
}

/*
** two functions for modifying an existing psys by adding
** monomials.
**  pysy_new_mon creates a new monomial in the
**      new_mon field, and sets curr_mon to point to it.
**      It can then be filled with data
**  once monomial has been filled with data psys_save_mon
**  will add it to a specified equation. (should not be used
**  while iterating through monomial list).
*/

int psys_init_mon(psys sys){
    if ((Mon_curr(sys)=(Sys_Mon_New(sys)))!=0) return TRUE;
    else return FALSE;
}

void psys_save_mon(psys sys, int i){
   int tmp=0,j,idx;
   /* calculate  degree of new monomial */
   for(j=1;j<=Sys_N(sys);j++) tmp+=Mon_exp(sys,Sys_Mon_New(sys),j);
  
   if (Eqn_start(sys,i)==0){ /* if first monomial set degree */
         Eqn_deg(sys,i)=tmp;
         Mon_homog(sys,Sys_Mon_New(sys))=0;
   }
   else if (Eqn_deg(sys,i)>=tmp){
        Mon_homog(sys,Sys_Mon_New(sys))=Eqn_deg(sys,i)-tmp;
   }
   else {
     Mon_homog(sys,Sys_Mon_New(sys))=0;
     idx=Eqn_start(sys,i);
     while(idx!=0){
       Mon_homog(sys,idx)+=tmp-Eqn_deg(sys,i);
       idx=Mon_next(sys,idx);
     }
     Eqn_deg(sys,i)=tmp;
   }  
   tmp=Mon_next(sys,Sys_Mon_New(sys));
   Mon_next(sys,Sys_Mon_New(sys))=Eqn_start(sys,i);
   Eqn_start(sys,i)=Sys_Mon_New(sys);
   Sys_Mon_New(sys)=tmp;
   Sys_M(sys)++; 
   Eqn_size(sys,i)++;
}

/*
** Display functions
*/
extern Pring Def_Ring;


psys psys_fprint(FILE *fout,psys sys){
 int i,pct=0,mct=0;
 Blk_curr(sys)=1;
  do {
    Eqn_curr(sys)=Blk_start(sys,Blk_curr(sys));
    do { 
      if (pct++==0)
#ifdef LOG_PRINT
 fprintf(fout,"< ")
#endif
;
      else 
#ifdef LOG_PRINT
fprintf(fout,",\n  ")
#endif
;
      Mon_curr(sys)=Eqn_start(sys,Eqn_curr(sys));
      mct=0;
      do {
        if (mct++!=0) 
#ifdef LOG_PRINT
fprintf(fout," + ")
#endif
;        if (Mon_coefR(sys,Mon_curr(sys))!=0||
             Mon_coefI(sys,Mon_curr(sys))!=0){
          if (Mon_coefI(sys,Mon_curr(sys)) != 0.0)
             fprintf(fout,"(");
          fprintf(fout,"%g", Mon_coefR(sys,Mon_curr(sys)));
          if (Mon_coefI(sys,Mon_curr(sys)) != 0.0) {
	    if (Mon_coefI(sys,Mon_curr(sys))>=0.0)fprintf(fout," + ");
	    fprintf(fout,"%g*I)", Mon_coefI(sys,Mon_curr(sys)));
	  }
          for(i=1;i<=Sys_N(sys);i++){
            if (Mon_exp(sys,Mon_curr(sys),i)!=0){
              fprintf(fout," %s^%d",ring_var(Def_Ring,i-1),
                               Mon_exp(sys,Mon_curr(sys),i));
            }
          }
          if (Mon_defv(sys,Mon_curr(sys))!=0){
             fprintf(fout," %s^%d",ring_def(Def_Ring),
                                     Mon_defv(sys,Mon_curr(sys)));
          }
          /* fprintf(fout,"\n"); */
        }
       }
       while((Mon_curr(sys)=Mon_next(sys,Mon_curr(sys)))!=0);
     }
     while(++Eqn_curr(sys)<Blk_start(sys,Blk_curr(sys)+1));
  }
  while((++Blk_curr(sys))<=Sys_R(sys));
 fprintf(fout," >\n");
return sys;
}   


/*
**  Accessors
*/

int psys_d(psys sys) {return sys->N;}


double *psys_coef_real(psys sys){
 return &(Mon_coefR(sys,Mon_curr(sys)));
}

double *psys_coef_imag(psys sys){
   return &(Mon_coefI(sys,Mon_curr(sys)));
}

int *psys_exp(psys sys,int d){ 
   return &(Mon_exp(sys,Mon_curr(sys),d));
 }

int *psys_homog(psys sys){
   return &((Mon_homog(sys,Mon_curr(sys))));
}

int *psys_def(psys sys){ 
   return &(Mon_defv(sys,Mon_curr(sys)));
}

/*
** Iteration Functions
*/

int psys_start_poly(psys sys){
      Eqn_curr(sys)=1;
      return TRUE;
}

int psys_next_poly(psys sys){ 
     if (++(Eqn_curr(sys))<=Sys_N(sys)) return TRUE;
        else return FALSE;
}

int psys_start_block(psys sys){
      Blk_curr(sys)=1;
      return TRUE;
}

int psys_next_block(psys sys){
     if (++Blk_curr(sys)<=Sys_R(sys)) return TRUE;
        else return FALSE;
}

int psys_set_eqno(psys sys, int eq){
  return (Eqn_curr(sys)=eq);
}

int psys_eqno(psys sys){
  return Eqn_curr(sys);
}

int psys_bkno(psys sys){
  return Blk_curr(sys);
}

int psys_Bstart_poly(psys sys,int i){         
      Eqn_curr(sys)=Blk_start(sys,i);       
      return TRUE;                                   
}                         
 
int psys_Bnext_poly(psys sys,int i){
     if (++(Eqn_curr(sys))<Blk_start(sys,i+1)) return TRUE;
        else return FALSE;                           
}                         

int psys_start_mon(psys sys){
    Mon_curr(sys)=Eqn_start(sys,Eqn_curr(sys));
    if (Mon_curr(sys)!=0) return TRUE;
    return FALSE;
}

int psys_next_mon(psys sys){
      if ((Mon_curr(sys)=Mon_next(sys,Mon_curr(sys)))!=0)return TRUE;
      else return FALSE;
 }

/*
** psys_type -- compute type vector for psys
**              (warning must have block structure set up)
*/
Imatrix psys_type(psys sys){
  Imatrix T;
  int i;
  T=Ivector_new(Sys_R(sys));
  for(i=1;i<=Sys_R(sys);i++) *IVref(T,i)=Blk_size(sys,i);
  return T;
}

/*
** psys_aux  find total number of monomials
*/
void **psys_aux(psys sys){
  return &(Mon_aux(sys,Mon_curr(sys)));
}

/*
** psys_size  find total number of monomials
*/
int psys_size(psys sys){
  return Sys_M(sys);
}

/* 
** finds number of consecutive equations starting with current
** equation share a common support. 
*/
int psys_r(psys sys){ return Sys_R(sys);}
int psys_block_size(psys sys){
    return Blk_size(sys,Blk_curr(sys));
}
/*
** finds number of monomials in current equation
*/
int psys_eq_size(psys sys){                 
    return Eqn_size(sys,Eqn_curr(sys));
} 

int *psys_block_start(psys sys,int i){return &(Blk_start(sys,i));}

/*
** copying a psys
*/
psys psys_copy(psys sys)
  { int i;
    psys res;
    res=psys_new(psys_d(sys),psys_size(sys),psys_r(sys));
    for(i=1;i<=psys_r(sys);i++){
       *psys_block_start(res,i)=*psys_block_start(sys,i);
    }
    FORALL_POLY(sys,
      FORALL_MONO(sys,
        psys_init_mon(res);
        *psys_aux(res)=*psys_aux(sys); 
        *psys_coef_real(res)=*psys_coef_real(sys);
        *psys_coef_imag(res)=*psys_coef_imag(sys);
        *psys_def(res)=*psys_def(sys);
        *psys_homog(res)=*psys_homog(sys);
        for(i=1;i<=psys_d(sys);i++) *psys_exp(res,i)=*psys_exp(sys,i);
        psys_save_mon(res,psys_eqno(sys));
      )
    )
    return res;
 }

/* end psys_def.c */

/**************************************************************************/
/****************** implementation code from psys_eval.c ******************/
/**************************************************************************/

/*
** psys_eval
**       input: psys sys (psys)
**              xpnt   X    (xpnt)
**             Dvector Y
**    effects: Y is resized (if nescessary) and filled with values
**             of polynomials of sys at X.
**    output: Y
**    error conditions: if X and sys are not comatable abort program
*/        
#define Y(i) (DVref(Y,i))
Dvector psys_eval(psys sys, xpnt X, Dvector Y){
  int n,i,eqno=0;
  fcomplex pval,tmp;
  n=psys_d(sys);
  if(xpnt_n(X)!=n) bad_error("bad argument in psys_eval");
  if(Y==0 || DVlength(Y)!=2*n) Y=Dmatrix_resize(Y,1,2*n);