pelpsys.cc

Gambit 是一个游戏库理论软件

  for(i=1;i<=2*n;i++) Y(i)=0;
  FORALL_POLY(sys,
    eqno++;
    pval=Complex(0.0,0.0);
    FORALL_MONO(sys,
      tmp=Complex(*psys_coef_real(sys),*psys_coef_imag(sys));
      for(i=1;i<=n;i++){
        tmp=Cmul(tmp,Cpow(xpnt_xi(X,i),*psys_exp(sys,i)));
      }
      tmp=Cmul(tmp,Cpow(xpnt_h(X),*psys_homog(sys)));
      tmp=RCmul(pow(xpnt_t(X),(long)(*psys_def(sys))),tmp);
      pval=Cadd(tmp,pval);
    )
    Y(2*eqno-1)=pval.r;
    Y(2*eqno)=pval.i;
  )
  return Y;
}
#undef Y

#define J(i,j) (DMref(J,i,j))
Dvector psys_jac(psys sys, xpnt X, Dmatrix J){
  int n,i,j,eqno=0;
  fcomplex tmp;
  n=psys_d(sys);
  if(xpnt_n(X)!=n) bad_error("bad argument in psys_eval");
  if(DMrows(J)!=2*n || DMcols(J)!=2*n+3) J=Dmatrix_resize(J,2*n,2*n+3);

  for(i=1;i<=2*n;i++)for(j=1;j<=2*n+3;j++) J(i,j)=0;
  
  FORALL_POLY(sys,
    eqno++;
    FORALL_MONO(sys,
      /* dM/DH */
      if (*psys_homog(sys)!=0){
        tmp=Complex(*psys_coef_real(sys),*psys_coef_imag(sys));
        for(i=1;i<=n;i++){
          tmp=Cmul(tmp,Cpow(xpnt_xi(X,i),*psys_exp(sys,i)));
        }
        tmp=RCmul(*psys_homog(sys),tmp);
        tmp=Cmul(tmp,Cpow(xpnt_h(X),*psys_homog(sys)-1));
        if (*psys_def(sys)!=0)tmp=RCmul(pow(xpnt_t(X),
					    (long)(*psys_def(sys))),tmp);
        J(2*eqno-1,1)+=tmp.r;
        J(2*eqno  ,2)+=tmp.r;
        J(2*eqno-1,2)-=tmp.i;
        J(2*eqno  ,1)+=tmp.i;
      }
     /* DM/DT */
      if (*psys_def(sys)!=0){
        tmp=Complex(*psys_coef_real(sys),*psys_coef_imag(sys));
        for(i=1;i<=n;i++){
          tmp=Cmul(tmp,Cpow(xpnt_xi(X,i),*psys_exp(sys,i)));
        }
        tmp=Cmul(tmp,Cpow(xpnt_h(X),*psys_homog(sys)));
        tmp=RCmul(pow(xpnt_t(X),(long)(*psys_def(sys)-1)),tmp);
        tmp=RCmul(*psys_def(sys),tmp);
        J(2*eqno-1,2*n+3)+=tmp.r;
        J(2*eqno  ,2*n+3)+=tmp.i;                          
      }
     /* DM/DXi*/
     for(j=1;j<=n;j++){
       if (*psys_exp(sys,j)!=0){
          tmp=Complex(*psys_coef_real(sys),*psys_coef_imag(sys));
          tmp=RCmul(*psys_exp(sys,j),tmp);
          tmp=Cmul(tmp,Cpow(xpnt_xi(X,j),(*psys_exp(sys,j))-1));
          for(i=1;i<=n;i++){
           if (i!=j) tmp=Cmul(tmp,Cpow(xpnt_xi(X,i),*psys_exp(sys,i)));
          }
          tmp=Cmul(tmp,Cpow(xpnt_h(X),*psys_homog(sys)));
          if (*psys_def(sys)!=0)tmp=RCmul(pow(xpnt_t(X),
					      (long)(*psys_def(sys))),tmp);
          J(2*eqno-1,2*j+1)+=tmp.r;
          J(2*eqno  ,2*j+2)+=tmp.r;
          J(2*eqno-1,2*j+2)-=tmp.i;
          J(2*eqno  ,2*j+1)+=tmp.i;
       }
     }
    )
  )
  return J;
}                                     
#undef J

/*
** psys_abs
**       input: psys sys (psys)
**              xpnt   X    (xpnt)
**             Dvector Y
**    effects: Y is resized (if nescessary) and filled with absolute
**             values of polynomials of sys at evaluated at X.
**    output: Y
**    error conditions: if X and sys are not comatable abort program
*/                

double psys_abs(psys sys, xpnt X){
  int n,i;
  double d;
  fcomplex pval,tmp;
  n=psys_d(sys);
  if(xpnt_n(X)!=n) bad_error("bad argument in psys_eval");
  d=0;
  FORALL_POLY(sys,
    pval=Complex(0.0,0.0);
    FORALL_MONO(sys,
      tmp=Complex(*psys_coef_real(sys),*psys_coef_imag(sys));
      for(i=1;i<=n;i++)
        tmp=Cmul(tmp,Cpow(xpnt_xi(X,i),*psys_exp(sys,i)));
      tmp=Cmul(tmp,Cpow(xpnt_h(X),*psys_homog(sys)));
      tmp=RCmul(pow(xpnt_t(X),(long)(*psys_def(sys))),tmp);
      pval=Cadd(tmp,pval);
    )
    d+=pval.r*pval.r+pval.i*pval.i;
  )
  return sqrt(d);
}

/*
** moment_sys
**      input: a pvector P
**             a xpnt X 
**                        
**    output:  a Dvector M representing a point on the newton
**             polytope of the system, the sum of the images of
**             the moment maps defined by the supporting Aset of
**             the system P.
**
**    note if M allready has enough space allocated to it will be
**         reused.
*/
#define M(i) DVref(M,i)
#define tmp(i) DVref(tmp,i)
Dvector psys_moment(psys sys,xpnt X, Dvector M){
   int i,n;
   double den,t;
   Dmatrix tmp;
   n=psys_d(sys);
   M=Dmatrix_resize(M,1,n);
   tmp=Dmatrix_new(1,n);
   for(i=1;i<=n;i++) M(i)=0.0;

   FORALL_BLOCK(sys,
     psys_Bstart_poly(sys,psys_bkno(sys));
     for(i=1;i<=n;i++) tmp(i)=0.0;
     den=1.0;
     FORALL_MONO(sys,
       t=1.0;
       for(i=1;i<=n;i++) t*=pow(Cabs(xpnt_xi(X,i)),(long)(*psys_exp(sys,i)));
       t*=pow(Cabs(xpnt_h(X)),(long)(*psys_homog(sys)));
       for(i=1;i<=n;i++) tmp(i)+=(*psys_exp(sys,i)*t);
       den+=t;
     );
     for(i=1;i<=n;i++) M(i)=tmp(i)/den;
   )
Dmatrix_free(tmp);
return M;}
#undef M
#undef tmp

/* end psys_eval.c */

/**************************************************************************/
/******************* implementation code from psys_hom.c ******************/
/**************************************************************************/

static Imatrix Norm=0;

node psys_hom(psys sys, node point_list, int tweak){
  node ptr=point_list;   

#ifdef PSYS_HOM_PRINT
  fprintf(stdout /* was psys_logfile */,"S starting continuation:\n");
  fprintf(stdout /* was psys_logfile */,"N");Imatrix_fprint(stdout /* was psys_logfile */,Norm);
  fprintf(stdout /* was psys_logfile */,"P"); psys_fprint(stdout /* was psys_logfile */,sys);
#endif

  if (Cont_Alg==USE_HOMPACK){
     init_hom(sys);
     while(ptr!=0){
       HPK_cont((Dmatrix)Car(Car(ptr)), tweak);
       ptr=Cdr(ptr);
     }
  }
  else {
    while(ptr!=0){

#ifdef PSYS_HOM_PRINT
      fprintf(stdout /* was psys_logfile */,"P"); psys_fprint(stdout /* was psys_logfile */,sys);
      fprintf(stdout /* was psys_logfile */,"X"); Dvector_fprint(stdout /* was psys_logfile */,(Dmatrix)Car(Car(ptr)));
#endif

      psys_cont((Dmatrix)Car(Car(ptr)),sys);
      ptr=Cdr(ptr);
    }
  }

#ifdef PSYS_HOM_PRINT
  fprintf(stdout /* was psys_logfile */,"Q ending continuation\n");
#endif

  return point_list;
}

node psys_solve(psys sys, Imatrix norm, int tweak){
  psys start,hom;
  node sols=0;
  Norm=norm;
  hom=psys_norm_sub(psys_copy(sys),norm);
  start=psys_lead(hom);
  sols=psys_hom(hom,psys_binsolve(start),tweak);
  Norm=0;
  psys_free(start);
  psys_free(hom);
  return sols;
}

/* end psys_hom.c */

/**************************************************************************/
/******************* implementation code from psys_ops.c ******************/
/**************************************************************************/

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

/*
** psys_lead(psys sys)
**      Input:  a psys.
**     Output:  (copies) of lowest terms of pys, with largest monomial
**              gcd removed.
**     side effects:
**     error conditions:
*/
psys psys_lead(psys sys)
  { int i,t,n,m;
    psys res;
    n=psys_d(sys);
    res=psys_new(n,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,
    /*
    ** find lowest exponent m of deformation parameter in current
    ** polynomial of input system sys - and copy monomials having
    ** this value over to the output system res
    */
      psys_start_mon(sys);
      m=(*psys_def(sys));
      if (psys_next_mon(sys)==TRUE){
        do{
           if (m>(t=*psys_def(sys))) m=t;
        }
        while(psys_next_mon(sys)==TRUE);
      }

      FORALL_MONO(sys,
        if (*psys_def(sys)==m){
           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);
           for(i=1;i<=n;i++) *psys_exp(res,i)=*psys_exp(sys,i);
           *psys_homog(res) = *psys_homog(sys);
           *psys_def(res)=0;
           psys_save_mon(res,psys_eqno(sys));
        }
      )
    )
    return psys_saturate(res);
  }

/*
** psys_saturate
**       input: psys 
**      output: psys with largest common monomial factor removed.
**      side effects: original psys is equal to output psys.
**      error contitions:
*/
#define V(i) (*IVref(V,i))
psys psys_saturate(psys sys)
  { int i,t,n;
    Ivector V;
    n=psys_d(sys);
    V=Ivector_new(n);
    FORALL_POLY(sys,
     /*
     ** find largest common monomial factor for all
     ** monomials in res (stored by vector V of exponents)
     ** WARNING: will also need to check homog variable
     */
      psys_start_mon(sys);
      for(i=1;i<=n;i++) V(i)=(*psys_exp(sys,i));
      if (psys_next_mon(sys)==TRUE){
        do{
          for(i=1;i<=n;i++){
            if ((t=*psys_exp(sys,i))<V(i)) V(i)=t;
          }
        }
        while(psys_next_mon(sys)==TRUE);
      }
      /*
      ** devide res by common monomial factor (with exponents V)
      */
      FORALL_MONO(sys,
           for(i=1;i<=n;i++) *psys_exp(sys,i)-=V(i);
      )
    )
    Ivector_free(V);
    return sys;
 }
#undef V

/*
** psys_lift:
**    input: a system S (psys)
**           a lifting value l (int)
**    output: original system with deformation parameter for each 
**            monomial set to l
**    side effects: original system equals output system.
*/
 psys psys_lift(psys sys, int l){
  FORALL_POLY(sys,FORALL_MONO(sys,*psys_def(sys)=l;));
 return sys;
 }

/*
** psys_norm_sub
**      input: a polynomial systm sys (psys)
**             a normal           norm (Ivector)
**      output: original system after transform
**              x_i->x_i*t^n_i. (and then removing common powers of t)
**      sidefects: original system equals output system
*/
#define norm(i) (*IVref(norm,i))
psys psys_norm_sub(psys sys,Ivector norm)
  { int i,t,n;
  int m = 0;
  int tog = 0;
    n=psys_d(sys);
    FORALL_POLY(sys,
      tog=0;
      FORALL_MONO(sys,
        t=(norm(n+1))*(*psys_def(sys));
        for(i=1;i<=psys_d(sys);i++)
            t+=(*psys_exp(sys,i))*norm(i);
        *psys_def(sys)=t;
        if (tog==0) { m=t; tog=1;}
        else if (m>t) m=t;
      )
      FORALL_MONO(sys,*psys_def(sys)-=m;)
    )
    return sys;
  }
#undef norm


/* end psys_ops.c */


/**************************************************************************/
/****************** implementation code from psys_scale.c *****************/
/**************************************************************************/

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

/*
** scale.c             Birk Huber                    created 4-8-1995
** All Rights Reserved
**
** Very simple program for polynomial system scaling.
**
**   Given a system F_1,...,F_n of polynomial equations
**   in variables X_1,...,X_n  calculates constants d_i
**   and c_i so that the new equations
**      10^{d_i}F_i(10^{c_1}X_1,...,10^{c_n}X_n)
**   will have coefitients as close to one as possible (i.e. the 
**   two norm of the vector of logarithms will be minimized).
**
**   The book on continuation by Morgan describes a similar algorithm,
**   in which variation among the coeficients is also minimized
**   explicitly. This code does not take speciall steps to also 
**   minimize the diferences -- It is sort of taken care of allready
**   by the least squares problem allready. (might add it later need 
**   only new rows in matrix).
*/ 

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

/* 
** Array and Polynomial data acces macroes
*/
#define LHS(i,j) (DMref(LHS,i,j))
#define RHS(i)   (DVref(RHS,i))
/* #define   X(i)   (DVref(X,i)) REDEFINITION */
#define CEXP(i)  (*psys_exp(Sys,i))
#define CCR      (*psys_coef_real(Sys))
#define CCI      (*psys_coef_imag(Sys))

/*
** Scale
**      Input: a polynomial system Sys (n-variables n-equations).
**     Output: Vector whoose ith coordinate is the inverse 10^{-c_i}
**             of the variable scaling factor for the i-th variable
**     Side Effect: Sys gets scaled.
*/

Dvector psys_scale(psys Sys){
 int i,j,k,n,row=1,m=0,eqno=0;
 Dmatrix LHS,RHS,X;
 double t,t1,t2,s,s1,s2;

 n=psys_d(Sys);
 m=psys_size(Sys);

/*
** Set up least squares problem: find X minimizing |(LHS)(X)-(RHS)|
**
**  Under scaling 10^d_i*F_i(x_1^e_1,...,x_n^e_n) each monomial m=a*x_1^e_1*...*x_n^e_n
**   becomes 10^{d_i}*a*10^{c_1*e_1+...+c_n*e_n}x_1^{e_1}*...x_n^{e_n}
**   taking logarithms the goal that the new coeficient should be equal to one 
**   becomes a linear condition  d_i+c_1*e_1+...c_n*e_n = -log(a)
**   which gets writen out as a matrix equation for the unknown
**         X=[d_1,...,d_n,c_1,...,c_n]. 
*/
 RHS=Dvector_new(m);
 LHS=Dmatrix_new(m,2*n);
 X=Dvector_new(2*n);

 eqno=0;
 psys_start_poly(Sys);
 do {
   eqno++;
   psys_start_mon(Sys);
   do {
     for(j=1;j<=n;j++) LHS(row,j)=0;
     LHS(row,eqno)=1.0;
     for(j=1;j<=n;j++) LHS(row,n+j)=CEXP(j);
     RHS(row)=-1.0*log10(sqrt(CCR*CCR+CCI*CCI)); 
     row++;
   }
   while(psys_next_mon(Sys)==TRUE);
 }
 while(psys_next_poly(Sys)==TRUE);

/* 
** Use Givens rotations to triangularize LHS,i.e.LHS becomes Q*LHS.(triangular) 
** and also apply same givens rotations to RHS  i.e. RHS becomes Q*RHS.
*/
 for (i=1;i<=2*n;i++){
   for (k=i+1;k<=m;k++){
     if (LHS(k,i)!=0.0){
         /* compute rotation */
         s=sqrt(LHS(i,i)*LHS(i,i)+LHS(k,i)*LHS(k,i));
         s1=LHS(i,i)/s;
         s2=LHS(k,i)/s;
         /* apply rotation to LHS */
         for(j=1;j<=2*n;j++) {
             t1=LHS(i,j);
             t2=LHS(k,j);
             LHS(i,j)=s1*t1+s2*t2;
             LHS(k,j)=-s2*t1+s1*t2;
         }
         /* apply rotation to RHS */
         t1=RHS(i);
         t2=RHS(k);
         RHS(i)=s1*t1+s2*t2;
         RHS(k)=-s2*t1+s1*t2;
     }
   }
 }

/* 
** Back solve   R*X=Y
**   R=top square portion of LHS. (2nx2nupper triangular)
**   Y=top 2n entrees of RHS.
**   results in X which solves original least squares problem.
*/
for(j=2*n;j>=1;j--){
  t=RHS(j);
  for(i=j+1;i<=2*n;i++) t-=LHS(j,i)*X(i);
  t=t/LHS(j,j);
  X(j)=t;
}

/* 
** Actually perform scaling   (X=[d_1,...,d_n,c_1,...,c_n])
**  adjust coeficients to those of new system (exponents remain unchanged)
**         10^{d_i}F_i(10^{c_1}x_1,...,10^{c_n}x_n)
*/
 eqno=0;
 psys_start_poly(Sys);
 do {
   eqno++;
   psys_start_mon(Sys);
   do {
   t=X(eqno);
   for(j=1;j<=n;j++) t+=(X(n+j)*CEXP(j));  
   CCR*=pow(10.0,t);
   CCI*=pow(10.0,t);
   }
   while(psys_next_mon(Sys)==TRUE);
 }
 while(psys_next_poly(Sys)==TRUE);


/* 
** clean-up and return vector defining inverse scaling factors 
** to be used to get solutions for original problem.
*/
Dmatrix_free(LHS);
Dmatrix_free(RHS);
RHS=Dvector_new(n);
for(i=1;i<=n;i++) RHS(i)=pow(10.0,X(n+i));
Dmatrix_free(X);
return RHS;
}


/*
** Undefine Array and Polynomial data access macroes
*/
#undef LHS
#undef RHS
#undef X
#undef CEXP
#undef CCR
#undef CCI


/* end psys_scale.c */