pelutils.cc

Gambit 是一个游戏库理论软件

double rand_double(int low, int high){
      return (drand48()*(high-low)+low);
}


/**************************************************************************/
/********************** implementations from Dlist.c **********************/
/**************************************************************************/

void Dlist_empty(node L){ Dlist_first(L)=0;}

/*
** invariants: Dnode_next(Dnode_prev(pos)) ==pos 
**             for all pointers to Dlist node entrees.
**             (NOTE: Dnode_next(L) := Dlist_first(L), 
**                    where L is list header) 
**
**             Dnode_prev(Dnode_next(pos))=pos
**             for all non-zero pointers to Dlist node entrees 
**             and also for list header
**             
*/

node Dlist_add(node L, node data){
  node tmp=0;
  LOCS(2);
  PUSH_LOC(data);
  PUSH_LOC(tmp);
  tmp=node_new();
  Dnode_link(tmp)=node_new();
  Dnode_data(tmp)=data;
  Dnode_prev(tmp)=L;
  Dnode_next(tmp)=Dlist_first(L);
  if (Dnode_next(tmp)!=0) Dnode_prev(Dnode_next(tmp))=tmp;
  Dlist_first(L)=tmp; 
  POP_LOCS();
  return tmp;
}


node Dlist_del(node L, node pos){
  if (Dnode_next(pos)!=0) Dnode_prev(Dnode_next(pos))=Dnode_prev(pos);
  if (Dnode_prev(pos)!=L) Dnode_next(Dnode_prev(pos))=Dnode_next(pos);
  else Dlist_first(L)=Dnode_next(pos);
  return (node)Dnode_data(pos);
}

node Dlist_data(node pos){ return (node)Dnode_data(pos);}

/* end Dlist.c */


/**************************************************************************/
/******************** implementations from Dmatrix.c **********************/
/**************************************************************************/

/*--------------------------------------------------------------- 
  vector/matrix type  a linear array of int, whith auxilary info.
       *) the number of elements that can be stored is in elt[0]
       *) the current number of rows is in elt[1]
       *) the current number of collumbs is in elt[2]
The actual data are then stored in row major order from elt[3] on
---------------------------------------------------------------*/
#ifndef DMATRIX_FAST
struct Dmatrix_t {
    int store;
    int nrows;
    int ncols;
    double *coords;
};
#endif

/*-------------------------------------------------------------
 vector access macroes (which ignore any rows except for first)
-------------------------------------------------------------*/
#define Vstore(V)  ((V->store))	/* maximum #elts available */
#define Vlength(V) ((V->nrows))	/* actual #elts stored      */
#define Vref1(V,i) (((V->coords)[i-1]))	   /*acces ith elt (starting at 1) */
#define Vref0(V,i)  (((V->coords)[i]))	/*acces ith elt (starting at 0) */
#define Vref(V,i)  Vref1(V,i)

/*------------------------------------------------------------
 matrix access macroes
-------------------------------------------------------------*/
#define Mstore(V)  ((V->store))	/* maximum #elts available */
#define MMrows(V)  ((V->store/V->ncols))	/* maximum #rows          */
#define Mrows(V) ((V->nrows))	/* number rows stored */
#define Mcols(V) ((V->ncols))	/* number cols stored */
#define Mref1(V,i,j)(((V->coords)[(i-1)*(V->ncols)+j-1]))
#define Mref0(V,i,j)(((V->coords)[i*(V->ncols)+j]))
#define Mref(V,i,j)  Mref1((V),i,j)

#ifndef DMATRIX_FAST
double DMstore(Dmatrix M)
{
    return Mstore(M);
}
double DMMrows(Dmatrix M)
{
    return MMrows(M);
}
double DMrows(Dmatrix M)
{
    return Mrows(M);
}
double DMcols(Dmatrix M)
{
    return Mcols(M);
}
double *DMelts(Dmatrix M)
{
    return M->coords;
}
double *DMref_P(Dmatrix M, int i, int j)
{
    return &(Mref0(M, i, j));
}

#endif 

/*
   **   Constructor/Destructors for Dmatrixes
   ** 
   ** Dmatrix Dmatrix_free(int r, int c); 
   **       New Dmatrix cabable of holding r rows, and c collumbs.
   ** Dmatrix Dmatrix_new(Dmatrix V);
 */

Dmatrix Dmatrix_new(int r, int c)
{
    Dmatrix V;
    /*    void *mem_malloc(int); */
    V = (Dmatrix) mem_malloc(sizeof(struct Dmatrix_t));
    if (!V)
	bad_error("allocation failure in Dmatrix_new()");
    V->coords = (double *) mem_malloc(r * c * sizeof(double));
    if (!V)
	bad_error("allocation failure 2 in Dmatrix_new()");
    Mstore(V) = r * c;
    Mrows(V) = r;
    Mcols(V) = c;
    return V;
}

void Dmatrix_free(Dmatrix V)
{   /* void mem_free(); */
    if (V != 0 && V->coords != 0)
	mem_free((char *) (V->coords));
    if (V != 0)
	mem_free((char *) (V));
}
#undef mem_malloc
#undef mem_free
/*
   ** Dmatrix_resize(R,r,c)
   **   if R has enough storage to hold an rxc matrix resets
   **   row and columb entrees of r to r and c. otherwise
   **   frees R and reallocates an rxc matrix
 */
Dmatrix Dmatrix_resize(Dmatrix R, int r, int c)
{

    if (R == 0 || Mstore(R) < (r * c)) {
	if (R != 0)
	    Dmatrix_free(R);
	R = Dmatrix_new(r, c);
    } else {
	Mrows(R) = r;
	Mcols(R) = c;
    }
    return R;
}


/*
   **  Dmatrix_print(M):  print an Dmatrix
   **    if M is null print <<>> and return fail.
   **    otherwise print matrix and return true.
 */
Dmatrix Dmatrix_fprint(FILE *fout, Dmatrix M)
{
    int i, j;

    if (M == 0) {
#ifdef LOG_PRINT
	fprintf(fout,"<<>>")
#endif
;
	return 0;
    }
#ifdef LOG_PRINT
    fprintf(fout,"<")
#endif
;
    for (i = 1; i <= Mrows(M); i++) {
#ifdef LOG_PRINT
	fprintf(fout,"<")
#endif
;
	for (j = 1; j <= Mcols(M); j++) {
#ifdef LOG_PRINT
	    fprintf(fout,"%8.4f", Mref(M, i, j))
#endif
;
	    if (j < Mcols(M))
#ifdef LOG_PRINT	
	fprintf(fout,", ")
#endif
;
	}
#ifdef LOG_PRINT
	fprintf(fout,">")
#endif
;
	if (i < Mrows(M))
#ifdef LOG_PRINT
	    fprintf(fout,"\n")
#endif
;
    }
#ifdef LOG_PRINT
    fprintf(fout,">")
#endif
;
    return M;
}

/*
   ** matrix_add(M1,M2,&M3) -- Add two Dmatrixes:
   **  if M1, and M2 are incompatable (or null) complain and return false.
   **  if *M3 has too little storage (or is null) free *M3 if nescesary
   **                                            and create new storage.
 */
Dmatrix Dmatrix_add(Dmatrix M1, Dmatrix M2, Dmatrix * M3)
{
    int i, j;
    Dmatrix R = *M3;

    if (M1 == 0 || M2 == 0 || Mrows(M1) != Mrows(M2) || Mcols(M1) != Mcols(M2)) {
#ifdef LOG_PRINT
	fprintf(stderr, "matrix_add: dimensions dont match\n")
#endif
;
	return 0;
    }
    Dmatrix_resize(R, Mrows(M1), Mcols(M1));
    for (i = 1; i <= Mrows(M1); i++)
	for (j = 1; j <= Mcols(M1); j++)
	    Mref(R, i, j) = Mref(M1, i, j) + Mref(M2, i, j);
    M3 = &R;
    return R;
}


/*
   **Dmatrix_mull(M1,M2,&M3) -- multiply two Dmatrixes:
   ** if M1, and M2 are incompatable (or null) complain and return false.
   ** if *M3 has too little storage (or is null) free *M3 if nescesary
   **                                           and create new storage.
   ** NOT USED YET
 */

Dmatrix Dmatrix_mull(Dmatrix M1, Dmatrix M2, Dmatrix * M3)
{
    int i, j, k;
    Dmatrix R = *M3;

    if (M1 == 0 || M2 == 0 || Mcols(M1) != Mrows(M2)) {
#ifdef LOG_PRINT
	fprintf(stderr, "Dmatrix_mull: incompatible matrices\n")
#endif
;
	return 0;
    }
    Dmatrix_resize(R, Mrows(M1), Mcols(M2));
    for (i = 1; i <= Mrows(M1); i++)
	for (j = 1; j <= Mcols(M2); j++) {
	    Mref(R, i, j) = 0;
	    for (k = 1; k <= Mcols(M1); k++)
		Mref(R, i, j) += Mref(M1, i, k) * Mref(M2, k, j);
	}
    M3 = &R;
    return R;
}

/*
   ** Dmatrix_dot(M1,M2,&M3) -- calculate M1*Transpose(M2):                
   ** if M1, and M2 are incompatable (or null) complain and return false.
   ** if *M3 has too little storage (or is null) free *M3 if nescesary 
   **                                           and create new storage.   
 */
Dmatrix Dmatrix_dot(Dmatrix M1, Dmatrix M2, Dmatrix M3)
{
    int i, j, k;

    if (M1 == 0 || M2 == 0 || Mcols(M1) != Mcols(M2)) {
#ifdef LOG_PRINT
	fprintf(stderr, "Dmatrix_dot: incompatible matrices\n")
#endif
;
	return 0;
    }
    M3 = Dmatrix_resize(M3, Mrows(M1), Mrows(M2));

    for (i = 1; i <= Mrows(M1); i++)
	for (j = 1; j <= Mrows(M2); j++) {
	    Mref(M3, i, j) = 0;
	    for (k = 1; k <= Mcols(M1); k++)
		Mref(M3, i, j) += Mref(M1, i, k) * Mref(M2, j, k);
	}
    return M3;
}


/*
** Dmatrix_equal(M1,M2)  
**     return true if matrices represented by M1 and M2 are equal.
**     false otherwise.                     NOT TESTED YET
*/
int Dmatrix_equal(Dmatrix M1, Dmatrix M2)
{
    int i, j;
    if (M1 == 0 && M2 == 0)
	return 1;
    if (M1==0||M2==0||(Mrows(M1)!=Mrows(M2))||(Mcols(M1)!=Mcols(M2)))
        return 0;
    for (i = 1; i <= Mrows(M1); i++)
	for (j = 1; j <= Mcols(M1); j++)
	    if (Mref(M1, i, j) != Mref(M2, i, j))
		return 0;

    return 1;
}

/*
** Dmatrix_GQR
**    Use givens QR on a sqaure matrix:
**        after completion Q should be orthogonal, R triangular
**        and Q*A = original matrix A
*/
void Dmatrix_GQR(Dmatrix Q,Dmatrix A)
        /* A->ch<=A->rh=Q->rh=Q->ch */
   { int i,j,k;
     int r,c;
     double s,s1,s2;
     double t1,t2;
     r=Mrows(A);
     c=Mcols(A);
       for(i=1;i<=r;i++){
              for(k=1;k<=r;k++)Mref(Q,i,k)=0.0;
              Mref(Q,i,i)=1.0;}

       for (i=1;i<=c;i++)
         for (k=i+1;k<=r;k++)
              /* performing givens rotations to zero A[k][i] */
             if (Mref(A,k,i)!=0){
               s=sqrt(Mref(A,i,i)*Mref(A,i,i)+
                      Mref(A,k,i)*Mref(A,k,i));
               s1=Mref(A,i,i)/s;
               s2=Mref(A,k,i)/s;
               for(j=1;j<=c;j++) {
                          t1=Mref(A,i,j);
                          t2=Mref(A,k,j);
                          Mref(A,i,j)=s1*t1+s2*t2;
                          Mref(A,k,j)=-s2*t1+s1*t2;}
            /* actually doing givens row rotations on transpose Q */
               for(j=1;j<=r;j++){
                    t1=Mref(Q,j,i);
                    t2=Mref(Q,j,k);
                    Mref(Q,j,i)=s1*t1+s2*t2;
                    Mref(Q,j,k)=-s2*t1+s1*t2;}
             }
         }

/* 
** Dmatrix_Solve, Solve a square matrix in place
*/
#define LHS(i,j) Mref(LHS,i,j)
#define RHS(i) Vref(RHS,i)
void Dmatrix_Solve(Dmatrix LHS, Dmatrix RHS,int n){
 int i,j,k,m=n;
 double s,s1,s2,t1,t2;

 /*
 ** 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<=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<=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=n;j>=1;j--){
   t1=RHS(j);
   for(i=j+1;i<=n;i++) t1-=LHS(j,i)*RHS(i);
   t1=t1/LHS(j,j);
   RHS(j)=t1;
 }
}
#undef LHS
#undef RHS

/* end Dmatrix.c */

/**************************************************************************/
/********************* implementations from Imatrix.c *********************/
/**************************************************************************/

#ifndef min
#define min(i,j) ((i) < (j) ? (i): (j))
#endif
#ifndef FALSE
#define FALSE 0
#endif
#ifndef TRUE
#define TRUE 1
#endif


#ifndef IMATRIX_FAST 
/*--------------------------------------------------------------- 
  vector/matrix type  a linear array of int, whith auxilary info.
       *) the number of elements that can be stored is in elt[0]
       *) the current number of rows is in elt[1]
       *) the current number of collumbs is in elt[2]
The actual data are then stored in row major order from elt[3] on
---------------------------------------------------------------*/
struct Imatrix_t {
    int store;	/* maximum number of helts reserved */
    int topc;	/* effective number of columbs */
    int topr;	/* effective number of rows */
    int ncols;	/* number of elts between first elts of rows */
    int *elts;	/* ptr to first elt of block of space */
};

/*
**  a couple of functions which are defined in files that
**  I don't want to include right now, when things stableize
**      bad_error -- used mainly to stop things when a function
**                   is passed bad input, as things stablizes
**                   less catastrophic reactions may be introduced.
**      mem_alloc,
**      mem_free -- call malloc and free with some extra book-keeping.
*/
#endif

/*
void bad_error(char *);
char *mem_malloc(size_t);
void mem_free(char *);
*/

/*-------------------------------------------------------------
 vector access macroes (which ignore any rows except for first)
-------------------------------------------------------------*/

#define ImatrixVstore(V)  (((V)->store))	/* maximum #elts available */
#define ImatrixVlength(V) (((V)->topc))	/* actual #elts stored      */
#define ImatrixVref1(V,i)  (((V)->elts[i-1]))/* acces ith elt (starting at 1)*/
#define ImatrixVref0(V,i)  (((V)->elts[i]))/* acces ith elt (starting at 0) */
#define ImatrixVref(V,i)  ImatrixVref1(V,i)

/*------------------------------------------------------------
 matrix access macroes
-------------------------------------------------------------*/
#define ImatrixMstore(V)  (((V)->store))	/* maximum #elts available */
#define ImatrixMMrows(V)  (((V)->store/(V)->ncols))   /* maximum #rows    */
#define ImatrixMrows(V) (((V)->topr))	/* number rows stored */
#define ImatrixMcols(V) (((V)->topc))	/* number cols stored */
#define ImatrixMNcols(V) (((V)->ncols))	/* number cols stored */
#define ImatrixMref1(V,i,j) ((((V))->elts[(i-1)*((V)->ncols)+j-1]))
#define ImatrixMref0(V,i,j) (((V)->elts[(i*(V)->ncols)+j]))
#define ImatrixMref(V,i,j)  ImatrixMref1(V,i,j)


#ifndef IMATRIX_FAST
int IMstore(Imatrix M)
{
    return ImatrixMstore(M);
}
int IMMrows(Imatrix M)
{
    return ImatrixMMrows(M);
}
int IMrows(Imatrix M)
{
    return ImatrixMrows(M);
}
int IMcols(Imatrix M)
{
    return ImatrixMcols(M);
}
int *IMref1(Imatrix M, int i, int j)
{
    return &(ImatrixMref1(M, i, j));
}
#endif



/*
   **   Constructor/Destructors for Imatrixes
   ** 
   ** Imatrix Imatrix_free(int r, int c); 
   **       New Imatrix cabable of holding r rows, and c collumbs.
   ** Imatrix Imatrix_new(Imatrix V);
 */

Imatrix Imatrix_new(int r, int c)
{
    Imatrix V;
    V = (Imatrix) mem_malloc(sizeof(struct Imatrix_t));
    if (!V)
	bad_error("allocation failure 1 in Imatrix_new()");