pelutils.cc

Gambit 是一个游戏库理论软件

    V->elts = (int *) mem_malloc(r * c * sizeof(int));
    if (!V)
	bad_error("allocation failure 2 in Imatrix_new()");
    V->store = r * c;
    V->topc = V->ncols = c;
    V->topr = r;
    return V;
}

void Imatrix_free(Imatrix V)
{
    if (V != 0) {
	if (V->elts != 0)
	    mem_free((char *) (V->elts));
	mem_free((char *) (V));
    }
}

/*
   ** Imatrix_resize(R,r,c)
   **   Reset R to hold an r,by  c matrix.
   **   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
   ** DOES NOT PRESERVE INDECIES OF EXISTING DATA
 */
Imatrix Imatrix_resize(Imatrix R, int r, int c)
{

    if (R == 0 || ImatrixMstore(R) < (r * c)) {
	if (R != 0)
	    Imatrix_free(R);
	R = Imatrix_new(r, c);
    } else {
	ImatrixMrows(R) = r;
	ImatrixMcols(R) = c;
	ImatrixMNcols(R) = c;
    }
    return R;
}


Imatrix Imatrix_submat(Imatrix R, int r, int c)
{

    if (R == 0 || c > ImatrixMNcols(R) || r > ImatrixMrows(R) * ImatrixMNcols(R)) {
	bad_error("bad subscripts or zero matrix in Imatrix_submat()");
    } else {
	ImatrixMrows(R) = r;
	ImatrixMcols(R) = c;
    }
    return R;
}


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

    if (M == 0) {
 #ifdef LOG_PRINT
       fprintf(fout,"<>\n")
#endif
;
        return 0;
    }
 #ifdef LOG_PRINT
   fprintf(fout,"<")
#endif
;
    for (i = 1; i <= ImatrixMrows(M); i++) {
        for (j = 1; j <= ImatrixMcols(M); j++) {
   #ifdef LOG_PRINT
         fprintf(fout,"%d", ImatrixMref(M, i, j))
#endif
;
            if (j < ImatrixMcols(M))  
#ifdef LOG_PRINT
                fprintf(fout,", ")
#endif
;
        }
        if (i < ImatrixMrows(M))
#ifdef LOG_PRINT
            fprintf(fout,";\n")
#endif
;
    }
#ifdef LOG_PRINT
    fprintf(fout,">\n")
#endif
;
    return M;
}

/*
** matrix_add(M1,M2,&M3) -- Add two Imatrixes:
** 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.
*/
Imatrix Imatrix_add(int i1,Imatrix M1,int i2, Imatrix M2, Imatrix M3)
{
    int i, j;
    Imatrix R = M3;

    if (M1 == 0 || M2 == 0 || ImatrixMrows(M1) != ImatrixMrows(M2) || ImatrixMcols(M1) != ImatrixMcols(M2)) {
#ifdef LOG_PRINT
	fprintf(stderr, "matrix_add: dimensions dont match\n")
#endif
;
	return 0;
    }
    R=Imatrix_resize(R, ImatrixMrows(M1), ImatrixMcols(M1));
    for (i = 1; i <= ImatrixMrows(M1); i++)
	for (j = 1; j <= ImatrixMcols(M1); j++)
	    ImatrixMref(R, i, j) =i1* ImatrixMref(M1, i, j) + i2*ImatrixMref(M2, i, j);
    return R;
}


/*
   **Imatrix_mull(M1,M2,&M3) -- multiply two Imatrixes:
   ** 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
 */

Imatrix Imatrix_mul(Imatrix M1, Imatrix M2, Imatrix M3)
{
    int i, j, k;
    Imatrix R;

    if (M1 == 0 || M2 == 0 || ImatrixMcols(M1) != ImatrixMrows(M2)) {
#ifdef LOG_PRINT	
fprintf(stderr, "Imatrix_mull: incompatible matrices\n")
#endif
;
	return 0;
    }
    if (M3!=M1&&M3!=M2) R=Imatrix_resize(M3, ImatrixMrows(M1), ImatrixMcols(M2));
      else R=Imatrix_new(ImatrixMrows(M1),ImatrixMcols(M2));
    for (i = 1; i <= ImatrixMrows(M1); i++)
	for (j = 1; j <= ImatrixMcols(M2); j++) {
	    ImatrixMref(R, i, j) = 0;
	    for (k = 1; k <= ImatrixMcols(M1); k++)
		ImatrixMref(R, i, j) += ImatrixMref(M1, i, k) * ImatrixMref(M2, k, j);
	}
    if (M3==M1&&M3==M2) Imatrix_free(M3);
    return R;
}

/*
   ** Imatrix_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.   
 */
Imatrix Imatrix_dot(Imatrix M1, Imatrix M2, Imatrix M3)
{
    int i, j, k;

    if (M1 == 0 || M2 == 0 || ImatrixMcols(M1) != ImatrixMcols(M2)) {
#ifdef LOG_PRINT
	fprintf(stderr, "Imatrix_dot: incompatible matrices\n")
#endif
;
	return 0;
    }
    M3 = Imatrix_resize(M3, ImatrixMrows(M1), ImatrixMrows(M2));

    for (i = 1; i <= ImatrixMrows(M1); i++)
	for (j = 1; j <= ImatrixMrows(M2); j++) {
	    ImatrixMref(M3, i, j) = 0;
	    for (k = 1; k <= ImatrixMcols(M1); k++)
		ImatrixMref(M3, i, j) += ImatrixMref(M1, i, k) * ImatrixMref(M2, j, k);
	}
    return M3;
}

/*
   ** dot_Ivector(M1,M2) -- vector dot product (integer valued)
   **  calculate dot product of first rows of M1, and M2.
   ** if M1, and M2 have are incompatable(or null) complain and return 0.
   ** NOT USED YET
 */
int dot_Ivector(Imatrix M1, Imatrix M2)
{
    int k, d = 0;
    if (M1 == 0 || M2 == 0 || ImatrixMcols(M1) != ImatrixMcols(M2)) {
#ifdef LOG_PRINT
	fprintf(stderr, "dot_Ivector: incompatible or null vectors\n")
#endif
;
	return 0;
    }
    for (k = 1; k <= ImatrixMcols(M1); k++)
	d += ImatrixVref(M1, k) * ImatrixVref(M2, k);
    return d;
}

/*
   ** Imatrix_equal(M1,M2)  
   **     return true if matrices represented by M1 and M2 are equal.
   **     false otherwise.  
 */
int Imatrix_equal(Imatrix M1, Imatrix M2)
{
    int i, j;
    if (M1 == 0 && M2 == 0)
	return TRUE;
    if (M1 == 0 || M2 == 0 || (ImatrixMrows(M1) != ImatrixMrows(M2)) || (ImatrixMcols(M1) != ImatrixMcols(M2))) {
	printf("in Imatrix_equal comparison failing\n");
	return FALSE;
    }
    for (i = 1; i <= ImatrixMrows(M1); i++)
	for (j = 1; j <= ImatrixMcols(M1); j++)
	    if (ImatrixMref(M1, i, j) != ImatrixMref(M2, i, j))
		return FALSE;

    return TRUE;
}


int Imatrix_order(Imatrix M1, Imatrix M2)
{
    int i, j, c = 0;
    if (M1 == 0 && M2 == 0)
	return TRUE;
    if (M1 == 0 || M2 == 0 || (ImatrixMrows(M1) != ImatrixMrows(M2)) || (ImatrixMcols(M1) != ImatrixMcols(M2))) {
	printf("in Imatrix_equal comparison failing\n");
	return FALSE;
    }
    for (i = 1; i <= ImatrixMrows(M1); i++)
	for (j = 1; j <= ImatrixMcols(M1); j++)
	    if ((c = ImatrixMref(M1, i, j) - ImatrixMref(M2, i, j)) != 0)
		return c;
    return 0;
}



int Imatrix_gcd_reduce(Imatrix M)
{
    int i, j, g = 0, gcd(int, int);

    if (M == 0)
	return 1;
    for (i = 1; i <= ImatrixMrows(M); i++)
	for (j = 1; j <= ImatrixMcols(M); j++) {
	    if (g == 0 && ImatrixMref(M, i, j) != 0)
		g = abs(ImatrixMref(M, i, j));
	    else if (ImatrixMref(M, i, j) != 0)
		g = abs(gcd(g, ImatrixMref(M, i, j)));
	}

    for (i = 1; i <= ImatrixMrows(M); i++)
	for (j = 1; j <= ImatrixMcols(M); j++)
	    ImatrixMref(M, i, j) = ImatrixMref(M, i, j) / g;

    return g;
}


int gcd(int a, int b)
{
    int c;
    if (b == 0)
	return a;
    c = a % b;
    if (c == 0)
	return b;
    else
	return gcd(b, c);
}

/*
** int Imatrix_rref(Imatrix N, int *det)
**   replaces N by its reduced row echelon form (upper triangular).
**   returns the rank of the  matrix, and puts the determinant
**   of the first full rank maximal minor in det.
**   ALGORITHM 9.1 of Algorithms for Computer Algebra,
**      Geddes,Czapor,Labahn
**  NOTE: I have actually changed this to reduce by the gcd when posible
**        to protect against over flow. (for the C-integer version).
*/
#define M(i,j) (ImatrixMref(N,i,j))

int Imatrix_rref(Imatrix N,int *det){
  int divisor=1,sign=1,r=1,p,k,j,i,g;

  for(k=1;k<=ImatrixMcols(N)&&r<=ImatrixMrows(N);k++){
    for(p=r; p<=ImatrixMrows(N)&&M(p,k)==0;p++); 
    if (p<=ImatrixMrows(N)) {
        for(j=k;j<=ImatrixMcols(N);j++){i=M(p,j); M(p,j)=M(r,j); M(r,j)=i;}
        if (r!=p) sign*=-1;
        for(i=r+1;i<=ImatrixMrows(N);i++){
          g=gcd(M(i,k),M(r,k)); 
          for(j=k+1;j<=ImatrixMcols(N);j++){
            M(i,j)=(M(r,k)*M(i,j)-M(r,j)*M(i,k))/divisor; 
          }
          M(i,k)=0;
        }
        divisor=M(r,k);
        r++;
    }
  }
  *det=sign*divisor;
  return r-1;
}
/*
** int Imatrix_backsolve(Imatrix N, Imatrix S)
**  finds a solution to a triangular system of equations with more
**  equtions then unknowns
*/
int Imatrix_backsolve(Imatrix N, Imatrix S){
 int p,j,g,k;


 if (N==0 || S==0 || IMcols(N)!=IMcols(S)||IMcols(N)<=IMrows(N)) 
     bad_error("bad dimensions in backsolve");
 
 /* 
 ** find first missing pivot (in row p;  p=IMrows(N)+1 if no missing pivots
 */
 for(p=1;p<=ImatrixMrows(N) && ImatrixMref(N,p,p)!=0;p++) /* null body*/;
 /* 
 ** now set up for a solution to the AX=0 problem with A a p-1xp matrix whoose
 ** first pxp minor is nonsingular (use simple minded back substitution)
 */
 for(j=1;j<=ImatrixMcols(N);j++) ImatrixVref(S,j)=0;
 ImatrixVref(S,p)=1;
 for (j=p-1;j>=1;j--){
   for(k=j+1;k<=p;k++) ImatrixVref(S,j)-=ImatrixMref(N,j,k)*ImatrixVref(S,k);
   g=gcd(ImatrixVref(S,j),ImatrixMref(N,j,j));
   ImatrixVref(S,j)=ImatrixVref(S,j)/g;
   g=(ImatrixMref(N,j,j)/g);
   if (ImatrixVref(S,j)!=0) for(k=j+1;k<=p;k++) ImatrixVref(S,k)*=g;
 }
return 1;
}
#undef M



/*
** Imatrix_hermite:
**    Input: S an n by m Imatrix.
**   Output: True 
**
**   Sidefects: U points to a unitary nxn matrix
**              S gets replaced by its hermite normal form (U*S)
**            
**  Method:
**    U = Inxn
**    c = 1
**    while (c<=n) do
**      find mc>=c such that abs(S(mc,c)) has minimum non-zero value
**      if (mc!=c) then  switch rows mc and c (in S, and U)
**      if (S(c,c)<0) then multiply row c by -1
**      if (S(c,c)!=0) then 
**        for i in c+1:n  add S(i,c)/S(c,c) times row c to row i;
**      end(if)
**      if S(i,c)==0 for all i in c+1 ... n set c=c+1
**    end(while)
**
*/

int Imatrix_hermite(Imatrix S, Imatrix U){
  int c=1,mc,i,j,m,n,done,sign;
  int t=0,mv=0;
  m=ImatrixMcols(S);
  n=ImatrixMrows(S);

  if (S==0 || U==0 || ImatrixMrows(U)!=n || ImatrixMrows(U)!=n)
     bad_error("Incompatable matrices in Imatrix_hermite");

  /* Initialize U to nxn identity */
  U=Imatrix_resize(U,n,n);
  for(i=1;i<=n;i++){
    for(j=1;j<=n;j++){
      if (i==j) ImatrixMref(U,i,j)=1;
      else ImatrixMref(U,i,j)=0;
    }     
  }

  while(c<=n){
    /* find minimum entry in col c */
    mv=abs(ImatrixMref(S,c,c));
    mc=c;
    for(i=c+1;i<=n;i++){
      t=abs(ImatrixMref(S,i,c)); 
      if(mv==0 || (mv>t && t!=0)){
           mv=t;
           mc=i;
      }
    }

    /* if nescesary pivot to put min in row c and multiply by+-1
       to ensure diagonal entry is positive */
    if (mc!=c||ImatrixMref(S,mc,c)<0){
      if (ImatrixMref(S,mc,c)<0) sign=-1;
      else sign=+1;
      for(j=c;j<=m;j++){
        t=ImatrixMref(S,mc,j);
        ImatrixMref(S,mc,j)=ImatrixMref(S,c,j);
        ImatrixMref(S,c,j)=sign*t;
      }
      for(j=1;j<=n;j++){
        t=ImatrixMref(U,mc,j);
        ImatrixMref(U,mc,j)=ImatrixMref(U,c,j);
        ImatrixMref(U,c,j)=sign*t;
      }
    }

    /* if collumb is not zero do a reduction step */
    done=TRUE;
    if (ImatrixMref(S,c,c)!=0){
      for(i=c+1;i<=n;i++){
        t=ImatrixMref(S,i,c)/ImatrixMref(S,c,c);
        for(j=c;j<=m;j++){
           ImatrixMref(S,i,j)-=t*ImatrixMref(S,c,j);
        }
        for(j=1;j<=n;j++){
           ImatrixMref(U,i,j)-=t*ImatrixMref(U,c,j);
        }
        if (ImatrixMref(S,i,c)!=0) done=FALSE;
      }
    }
    /* if all entrees of col c bellow row c are zero go tonext col */
    if (done==TRUE) c++;
  }
  return TRUE;
 }

Imatrix Imatrix_dup(Imatrix M,Imatrix N){
  Imatrix R;
  int i,j;
  if (M==0) return 0;
  if (M==N) bad_error("Imatrix_dup: source and destination equal");
  R=Imatrix_resize(N,ImatrixMrows(M),ImatrixMcols(M));
  for(i=1;i<=ImatrixMrows(M);i++)
     for(j=1;j<=ImatrixMcols(M);j++) ImatrixMref(R,i,j)=ImatrixMref(M,i,j);
  return R;
}


int Imatrix_is_zero(Imatrix M){
  int i,j;
  if (M==0) return TRUE;
  for(i=1;i<=IMrows(M);i++){
       for(j=1;j<=IMcols(M);j++){
         if (*IMref(M,i,j)!=0) return FALSE;
       }
  }
  return TRUE;
}

/* end Imatrix.c */

/**************************************************************************/
/********************** implementations from Lists.c **********************/
/**************************************************************************/

void bad_error(char *);
void mem_free(void *);

node list_push(node item,node *stack){
   *stack=Cons(item,*stack);
   return item;
} 
  
int list_length(node L){
  int ct=0;
  while(L!=0){
    ct++;
    L=Cdr(L);
  }
  return ct;
}

node list_pop(node *stack){
   node temp;
   temp=Car(*stack);
   *stack=Cdr(*stack);
   return temp;
}
  
node list_first(node list){
  return Car(list);
}

node list_rest(node list){
  return Cdr(list);
}

int list_empty(node list){
  if (list==0) return TRUE;
  else return FALSE;
}

/* 
** Inserts a node g onto a list L in order (according to a comparison
** function comp) -- IF There is allready and equivalent object on
** the list and uniq is TRUE it is not inserted, and a value of zero 
** is returned, otherwise it is inserted and a value of 1 is returned.
*/

int list_insert(node g, node * L, int (*comp) (node, node), int uniq)
{
  int flg = 1;
  node ptr = *L;
  LOCS(2);
  PUSH_LOC(*L);
  PUSH_LOC(g);

  if ((*L == 0) || ((flg = comp(g, Car(*L))) >= 0)) {
    if (uniq==FALSE || flg != 0) *L = Cons(g, *L);
  } 
  else {
    while ((Cdr(ptr) != 0) &&
      ((flg = comp(g, Car(Cdr(ptr))))< 0)) {
       ptr = (node) Cdr(ptr);
    }
    if (uniq==FALSE || flg != 0)
      node_set_ptr(ptr,Cons(g, Cdr(ptr)), NODE, RIGHT);
  }
  POP_LOCS();
  return flg;
}

/*