Game analysis program

Last updated on June 3, 2002.

The following program can be used to analyze two player constant sum games. It handles games with multiple states and with some types of recursion. I also have a downloadable version that does not have specials (<, > and &) escaped.

/*
 * Game Strategy Solver
 *
 * Written by Bruno Wolff III
 * December 10, 1989
 *
 * Changed "=-" to "= -" 12/26/91
 *
 * Made test for suitable pivots slightly larger than 0 rather than against
 * exactly 0. 6/4/96
 *
 * This program will for solve for the value of each state of a game
 * and for the optimal strategy for each state. Each state may have
 * its own set of strategies and state transitions.
 *
 * The input to the program is:
 *     The type of recursion allowed.
 *         (1 for no recursion, 2 for variable recursion, 3 for total recursion)
 *     The input format (1 for dense, 2 for sparse).
 *     The scale value to divide input values by.
 *     The number of states.
 *     If recursion allowed:
 *         The error limit to be used to check for completion.
 *         The maximum number of iterations.
 *     For each state:
 *         The number of row and number of column strategies.
 *         The immediate payoff table from the row point of view.
 *         The transition table to allowed states.
 *
 * The output from the program is:
 *    If recursion allowed:
 *        The number of iterations used.
 *        The maximum change in a state value in the last iteration.
 *    The value of each state.
 *    For total recursion:
 *        The value per turn of the game.
 *        The probability of being in each state.
 *    For each state:
 *        The optimal row strategy mix.
 *        The optimal column strategy mix.
 *
 * The problem is solved by treating each state as a seperate game
 * using the previous state values when accounting for state transitions.
 * Each game is then interpretted as a linear programming problem and
 * solved using the simplex method. The new state values are checked against
 * the old state values and the number of iterations is checked against the
 * maximum number of iterations to see if any more iterations are needed.
 * If more iterations are needed then new state values are calculated
 * using a multinomial Newton's method (unless I-F' is singular, in which
 * case the current new values are used) using the values of the transition
 * matrices (evaluted using the latest strategy mixes) as the entries of F'.
 * The new state values are then used as the old state values in the next
 * iteration. For total recursion an extra condition is added to I-F' so
 * that it will not be singular. The condition is that the value of each
 * state times the probability of being in that state sum to zero. In this
 * mode only one infinite loop of states is allowed.
 */

/*
 * Include files.
 */

#include <stdio.h>

/*
 * Declare function types.
 */

void main();
void gamsol();
int lineqs();
void badinp();

/*
 * Main program.
 */

void main()

{
int rec; /* type of recursion allowed */
int inp; /* format of input */
float scale; /* scale factor for input */
int nstate; /* number of states */
int miter; /* maximum number of iterations */
int iter; /* number of iterations */
float merr; /* maximum error */
float err; /* largest state value change in last iteration */
int *m; /* table of number of row startegies */
int *n; /* table of number of column strategies */
struct entry { /* sparse table entry */
  int estate; /* transition to state, 0 for payoff */
  float eval; /* value of entry */
  struct entry *enext; /* next entry in this state, row, column */
};
float ****dp; /* [state from][state to][row][col] dense payoff table */
struct entry ****sp; /* [state from][row][col] sparse payoff table */
struct entry *ent; /* scratch */
float *old; /* table of old state values */
float *val; /* table of current state values */
float **row; /* [state][strategy] row strategy mixes */
float **col; /* [state][strategy] column strategy mixes */
float **a; /* [row][col] scratch matrix used for gamsol and lineqs calls */
float *b; /* scratch vector used for lineqs call */
float *x; /* scratch vector used for lineqs call */
float v; /* scratch */
int kk; /* used to keep track of allowed state transitions */
int i,j,k,ii,jj; /* scratch */

  fprintf(stderr,"\nWhat type of recursion is allowed (1 none, 2 variable, 3 total)?\n");
  if (scanf("%d",&rec)<1) badinp();
  if (rec<1||rec>3) {
    fprintf(stderr,"\nThe type of recursion allowed must be 1, 2 or 3.\n");
    exit();
  }
  fprintf(stderr,"\nWhat input format will be used (1 dense, 2 sparse)?\n");
  if (scanf("%d",&inp)<1) badinp();
  if (inp<1||inp>2) {
    fprintf(stderr,"\nThe input format must be 1 or 2.\n");
    exit();
  }
  fprintf(stderr,"\nWhat is the scale factor for the input?\n");
  if (scanf("%f",&scale)<1) badinp();
  if (scale==0.) {
    fprintf(stderr,"\nThe scale factor must be nonzero.\n");
    exit();
  }
  fprintf(stderr,"\nHow many states does the game have?\n");
  if (scanf("%d",&nstate)<1) badinp();
  if (nstate < 1) {
    fprintf(stderr,"\nThe number of states can not be less than one.\n");
    exit();
  }
  kk=nstate; /* Assume any state transitons are allowed, change later if not */
  if (rec==1) {
    miter=1; /* Only one iteration needed if there is no recursion */
    merr=0.; /* This doesn't matter if there is no recursion */
  } else {
    fprintf(stderr,"\nWhat is the maximum error?\n");
    if (scanf("%f",&merr)<1) badinp();
    if (merr<=0.) {
      fprintf(stderr,"\nThe maximum error must be greater than zero.\n");
      exit();
    }
    fprintf(stderr,"\nWhat is the maximum number of iterations?\n");
    if (scanf("%d",&miter)<1) badinp();
    if (miter<1) {
      fprintf(stderr,"\nThe maximum number of iterations must be at least one.\n");
      exit();
    }
  }

/*
 * Allocate some space for tables.
 */

  m=(int *)malloc(sizeof(*m)*nstate); /* number of row strategies table */
  n=(int *)malloc(sizeof(*n)*nstate); /* number of column strategies table */
  old=(float *)malloc(sizeof(*old)*(nstate+2)); /* old value table */
  val=(float *)malloc(sizeof(*val)*(nstate+2)); /* new value table */
  row=(float **)malloc(sizeof(*row)*nstate); /* row strategy mixes */
  col=(float **)malloc(sizeof(*col)*nstate); /* column strategy mixes */
  b=(float *)malloc(sizeof(*b)*(nstate+1));
  x=(float *)malloc(sizeof(*x)*(nstate+1));
  if (inp==1) { /* Dense input space allocation */
    dp=(float ****)malloc(sizeof(*dp)*nstate);
    for (i=0; i<nstate; i++) {
      if (rec==1) kk=i; /* In case not all transitions are allowed */
      dp[i]=(float ***)malloc(sizeof(**dp)*(kk+1));
    }
  } else { /* Sparse input allocation */
    sp=(struct entry ****)malloc(sizeof(*sp)*nstate);
  }

/*
 * Get number of row and column strategies for each state.
 */

  for (i=0; i<nstate; i++) {
    fprintf(stderr,"\nWhat is the number of row and the number of column strategies for state %d?\n",i+1);
    if (scanf("%d%d",m+i,n+i)<2) badinp();
    if (m[i]<1) {
      fprintf(stderr,"\nThe number of row strategies must be at least one.\n");
      exit();
    }
    if (n[i]<1) {
      fprintf(stderr,"\nThe number of column strategies must be at least one.\n");
      exit();
    }
    row[i]=(float *)malloc(sizeof(**row)*m[i]); /* Space for row strategies */
    col[i]=(float *)malloc(sizeof(**col)*n[i]); /* Space for column strategies */
  }

/*
 * Enter tables using dense input format.
 */

  if (inp==1) {
    for (i=0; i<nstate; i++) {
      if (rec==1) kk=i; /* In case not all transitions are allowed */
      for (ii=0; ii<=kk; ii++) { /* Allocate space for each state */
	dp[i][ii]=(float **)malloc(sizeof(***dp)*m[i]);
	for (k=0; k<m[i]; k++) {
	  dp[i][ii][k]=(float *)malloc(sizeof(****dp)*n[i]);
	}
      }
      fprintf(stderr,"\nEnter %d x %d table of payoffs to the row player for state %d.\n",m[i],n[i],i+1);
      for (k=0; k<m[i]; k++) {
        for (j=0; j<n[i]; j++) {
          if (scanf("%f",dp[i][0][k]+j)<1) badinp();
	  dp[i][0][k][j]/=scale;
        }
      }
      for (ii=1; ii<=kk; ii++) {
        fprintf(stderr,"\nEnter %d x %d transition table from state %d to state %d.\n",m[i],n[i],i+1,ii);
        for (k=0; k<m[i]; k++) {
          for (j=0; j<n[i]; j++) {
            if (scanf("%f",dp[i][ii][k]+j)<1) badinp();
	    dp[i][ii][k][j]/=scale;
	    if (rec==3 && dp[i][ii][k][j]<0.) {
	      fprintf(stderr,"\nNegative transitions not allowed for total recursion.\n");
	      exit();
	    }
          }
        }
      }

/*
 * For total recursion make sure recursion adds up to one in every case.
 */

      if (rec==3) {
	for (k=0; k<m[i]; k++) {
	  for (j=0; j<n[i]; j++) {
	    v=0.;
	    for (ii=1; ii<=kk; ii++) v+=dp[i][ii][k][j];
	    if (v-1.<-1e-6 || v-1.>1e-6) {
	      fprintf(stderr,"\nRecursion of %f for state %d, row %d, column %d.\n",v,i+1,k+1,j+1);
	      exit();
	    }
	  }
	}
      }
    }
  }

/*
 * Enter tables using sparse format.
 */

  else {
    for (i=0; i<nstate; i++) { /* allocate space for each state */
      sp[i] = (struct entry ***)malloc(sizeof(**sp)*m[i]);
      for (k=0; k<m[i]; k++) {
	sp[i][k] = (struct entry **)malloc(sizeof(***sp)*n[i]);
	for (j=0; j<n[i]; j++) sp[i][k][j] = NULL;
      }
    }
    fprintf(stderr,"\nEnter payoff and transition table entries in the following formats:\n");
    fprintf(stderr,"State       Row  Column  0         Payoff to row player\n");
    fprintf(stderr,"State from  Row  Column  State to  Transition value\n");
    fprintf(stderr,"0 (this signifies the end of table entries)\n");
    while(1) {
      if (scanf("%d",&i)<1) badinp();
      if (i==0) break;
      if (i<0 || i>nstate) {
	fprintf(stderr,"\nBad state %d used for state from.\n",i);
	exit();
      }
      if (scanf("%d%d%d%f",&k,&j,&ii,&v)<4) badinp();
      if (k<1 || k>m[i-1]) {
	fprintf(stderr,"\nBad row %d used for state %d.\n",k,i);
	exit();
      }
      if (j<1 || j>n[i-1]) {
	fprintf(stderr,"\nBad column %d used for state %d.\n",j,i);
	exit();
      }
      if (ii<0 || ii>nstate || (rec==1 && ii>=i)) {
	fprintf(stderr,"\nBad state %d used for state to.\n",ii);
	exit();
      }
      if (v<0 && rec==3 && ii>0) {
	fprintf(stderr,"\nNegative transition in total recursion mode.\n");
	exit();
      }
      ent=sp[i-1][k-1][j-1];
      while (ent) {
	if (ent->estate==ii) {
	  fprintf(stderr,"\nDuplicate entry for state to %d, row %d, column %d, state from %d.\n",i,j,k,ii);
	  exit();
	}
	ent = ent->enext;
      }
      ent = (struct entry *)malloc(sizeof(*ent));
      ent->enext=sp[i-1][k-1][j-1];
      ent->eval=v/scale;
      ent->estate=ii;
      sp[i-1][k-1][j-1]=ent;
    }

/*
 * Make sure recursion adds up to one for total recursion.
 */

    if (rec==3) {
      for (i=0; i<nstate; i++) {
	for (k=0; k<m[i]; k++) {
	  for (j=0; j<n[i]; j++) {
	    v=0.;
	    ent=sp[i][k][j];
	    while (ent) {
	      if (ent->estate!=0) v+=ent->eval;
	      ent=ent->enext;
	    }
	    if (v-1.<-1e6 || v-1.>1e-6) {
	      fprintf(stderr,"\nRecursion of %f for state %d, row %d, column %d.\n",v,i+1,k+1,j+1);
	      exit();
	    }
	  }
	}
      }
    }
  }

/*
 * Allocate space for scratch matrix
 */

  ii=nstate+1;
  for (i=0; i<nstate; i++) {
    if (m[i]>=ii) ii=m[i]+1;
  }
  a=(float **)malloc(sizeof(*a)*ii);
  for (i=0; i<ii; i++) {
    j=nstate+1;
    for (k=0; k<nstate; k++) {
      if (m[k]>=i && n[k]>=j) j=n[k]+1;
    }
    a[i]=(float *)malloc(sizeof(**a)*j);
  }

/*
 * Use zero for initial state values.
 */

  old[0]=1.; /* value for payoff is special */
  for (i=1; i<=nstate+1; i++) old[i]=0.;

/*
 * Solve the game
 */

  for (iter=1; ; iter++) { /* keep going until the answer is good enough */

/*
 * Solve each state seperately.
 */

    for (i=0; i<nstate; i++) {
      if (rec==1) kk=i; /* recursion allowed? */

/*
 * Build parameters for calling gamsol.
 */

      if (inp==1) { /* dense input */
	for (j=0; j<n[i]; j++) {
	  for (k=0; k<m[i]; k++) {
	    v=dp[i][0][k][j];
	    for (ii=1; ii<=kk; ii++) v+=old[ii]*dp[i][ii][k][j];
	    a[k][j]=v;
	  }
	}
      }
      else { /* sparse input */
	for (j=0; j<n[i]; j++) {
	  for (k=0; k<m[i]; k++) {
	    v=0.;
	    ent=sp[i][k][j];
	    while(ent) {
	      v+=ent->eval*old[ent->estate];
	      ent=ent->enext;
	    }
	    a[k][j]=v;
	  }
	}
      }
      gamsol(m[i],n[i],a,row[i],col[i],val+i+1);
      if (rec==1) old[i+1]=val[i+1]; /* if no recursion value can be used immediately */
    }

/*
 * Calculate state probabilities for total recursion.
 * [ T    ]     [ ]
 * [F'-I 1] X = [0] solved for X gives the state probabilities.
 * [1    0]     [1]
 */

    if (rec==3) {
      if (inp==1) { /* dense input */
        for (i=0; i<nstate; i++) {
          for (ii=0; ii<nstate; ii++) {
            v=i==ii?-1.:0.;
            for (j=0; j<n[i]; j++) {
              for (k=0; k<m[i]; k++) {
                v+=row[i][k]*col[i][j]*dp[i][ii+1][k][j];
              }
            }
            a[ii][i]=v;
          }
        }
      }
      else { /* sparse input */
	for (i=0; i<nstate; i++) {
	  for (ii=0; ii<nstate; ii++) a[ii][i]=i==ii?-1.:0.;
	  for (j=0; j<n[i]; j++) {
	    for (k=0; k<m[i]; k++) {
	      ent=sp[i][k][j];
	      while (ent) {
		if (ent->estate!=0) a[ent->estate-1][i]+=row[i][k]*col[i][j]*ent->eval;
		ent=ent->enext;
	      }
	    }
	  }
	}
      }
      for (i=0; i<nstate; i++) {
	b[i]=0;
	a[i][nstate]=1.;
	a[nstate][i]=1.;
      }
      b[nstate]=1.;
      a[nstate][nstate]=0.;
      i=lineqs(nstate+1,a,b,x);
      if (i!=nstate+1) {
	fprintf(stderr,"\nSingular matrix for probability calculation in iteration %d.\n",iter);
	for (i=0; i<nstate; i++) x[i]=1./nstate;
      }
      val[nstate+1]=0.;
      for (i=0; i<nstate; i++) val[nstate+1]+=x[i]*val[i+1];
      for (i=1; i<=nstate; i++) val[i]-=old[nstate+1];
    }

/*
 * Check to see if we have done enough iterations.
 */

    err=0.;
    for (i=1; i<=nstate+1; i++) {
      v=val[i]-old[i];
      if (v<0.) v= -v;
      if (v>err) err=v;
    }
    if (rec==1||iter>=miter||err<=merr) break;

/*
 * Solve (I-F') X = Old - New, then set Old = Old - X. This will greatly speed
 * up convergence if (I-F') is not singular. If it is then set Old = New.
 */

    if (rec==2) { /* general recursion */
      if (inp==1) { /* dense input */
        for (i=0; i<nstate; i++) {
          b[i]=old[i+1]-val[i+1];
          for (ii=0; ii<nstate; ii++) {
            v=i==ii?1.:0.;
            for (j=0; j<n[i]; j++) {
              for (k=0; k<m[i]; k++) {
                v-=row[i][k]*col[i][j]*dp[i][ii+1][k][j];
              }
            }
            a[i][ii]=v;
          }
        }
      }
      else { /* sparse input */
	for (i=0; i<nstate; i++) {
	  b[i]=old[i+1]-val[i+1];
	  for (ii=0; ii<nstate; ii++) a[i][ii]=i==ii?1.:0.;
	  for (j=0; j<n[i]; j++) {
	    for (k=0; k<m[i]; k++) {
	      ent=sp[i][k][j];
	      while (ent) {
		if (ent->estate!=0) a[i][ent->estate-1]-=row[i][k]*col[i][j]*ent->eval;
		ent=ent->enext;
	      }
	    }
	  }
	}
      }
      i=lineqs(nstate,a,b,x);
      if (i==nstate) for (i=0; i<nstate; i++) old[i+1]=old[i+1]-x[i];
      else {
	fprintf(stderr,"\nSingular matrix in iteration %d.\n",iter);
	for (i=1; i<=nstate; i++) old[i]=val[i];
      }
    }

/*
 *       (    [F -1])
 * Solve (I - [P  0]) X = Old - New and then set Old = Old - X.
 * This will greatly speed up convergence if the matrix is not singular.
 * If it is set Old = New.
 */
    else { /* total recursion */
      if (inp==1) { /* dense input */
        for (i=0; i<nstate; i++) {
          for (ii=0; ii<nstate; ii++) {
            v=i==ii?1.:0.;
            for (j=0; j<n[i]; j++) {
              for (k=0; k<m[i]; k++) {
                v-=row[i][k]*col[i][j]*dp[i][ii+1][k][j];
              }
            }
            a[i][ii]=v;
          }
        }
      }
      else { /* sparse input */
	for (i=0; i<nstate; i++) {
	  for (ii=0; ii<nstate; ii++) a[i][ii]=i==ii?1.:0.;
	  for (j=0; j<n[i]; j++) {
	    for (k=0; k<m[i]; k++) {
	      ent=sp[i][k][j];
	      while (ent) {
		if (ent->estate!=0) a[i][ent->estate-1]-=row[i][k]*col[i][j]*ent->eval;
		ent=ent->enext;
	      }
	    }
	  }
	}
      }
      for (i=0; i<=nstate; i++) b[i]=old[i+1]-val[i+1];
      for (i=0; i<nstate; i++) {
	a[i][nstate]=1.;
	a[nstate][i]= -x[i]; /* copy state probabilities */
      }
      a[nstate][nstate]=1.;
      i=lineqs(nstate+1,a,b,x);
      if (i==nstate+1) for (i=0; i<=nstate; i++) old[i+1]=old[i+1]-x[i];
      else {
	fprintf(stderr,"\nSingular matrix in iteration %d.\n",iter);
	for (i=1; i<=nstate; i++) old[i]=val[i]+old[nstate+1]-val[nstate+1];
	old[nstate+1]=val[nstate+1];
      }
    }
  }

/*
 * Print out the results.
 */

  if (rec!=1) {
    printf("\nTermination after %d iterations.\n",iter);
    printf("The largest value change in the last iteration was %e\n",err);
  }
  if (rec==3) {
    printf("\nThe value per turn of the game is %.8f\n",val[nstate+1]);
    printf("\nThe game state probabilities and values are:\n");
    for (i=1; i<=nstate; i++) printf("State %4d    %.8f    %.8f\n",i,x[i-1],val[i]);
  }
  else {
    printf("\nThe game state values are:\n");
    for (i=1; i<=nstate; i++) printf("State %4d    %.8f\n",i,val[i]);
  }
  for (i=0; i<nstate; i++) {
    printf("\nThe active row strategies for state %d are:\n",i+1);
    for (k=0; k<m[i]; k++) {
      if (row[i][k]>=1.e-6) printf("Strategy %4d    %.8f\n",k+1,row[i][k]);
    }
    printf("\nThe active column strategies for state %d are:\n",i+1);
    for (j=0; j<n[i]; j++) {
      if (col[i][j]>=1.e-6) printf("Strategy %4d    %.8f\n",j+1,col[i][j]);
    }
  }
  printf("\n");
}

/*
 * Something about the input is messed up.
 */

void badinp()

{
  fprintf(stderr,"\nBad input format or end of file reached.\n");
  exit();
}

/*
 * Subroutine to solve a zero sum game using the simplex method.
 *
 * Written by Bruno Wolff III
 * June 13, 1989
 *
 * The game A is solved by minimizing the sum of X, with AX>=1, X>=0.
 * For this to be solvable the value of the game must be positive. So
 * that game A' is actually solved where A'=A-min(A)+1. The final answer
 * is given by X=X/(sum of X). V=1/(sum of X)+min(A)-1. The dual solution
 * is calculated at the same time and contains the other player's optimal
 * strategy.
 *
 * Parameters:
 * m      Number of row strategies (input)
 * n      Number of column strategies (input)
 * a      Payoff matrix (to row player) with an extra row and column (input/scratch)
 * r      m optimal row strategies (output)
 * c      n optimal column strategies (output)
 * v      Value of the game to the row player (output)
 */

void gamsol(m,n,a,r,c,v)

int m;
int n;
float **a; /* a has an extra row and column to be used for work space */
float *r;
float *c;
float *v; /* pointer so as to allow a value to be returned */

{
float s; /* minimum value of the a array minus one */
float p; /* best value of the pivot criterion so far */
float pc; /* curent pivot criterion */
float cp; /* best value of the pivot criterion in the current column */
int ip,jp; /* location of best pivot */
int icp; /* location of the best pivot in the current column */
int i,j; /* scratch */
int *rin; /* active row strategies */
int *rout; /* inactive row strategies */
int *cin; /* active column strategies */
int *cout; /* inactive column strategies */

/*
 * Check reasonableness of problem size.
 */

  if (m<1) {
    fprintf(stderr,"\ngamsol error - m must be at least 1.\n");
    exit();
  }
  if (n<1) {
    fprintf(stderr,"\ngamsol error - n must be at least 1.\n");
    exit();
  }

/*
 * Make all payoffs positive. This assures the value of the game is
 * positive.
 */

  s=a[0][0]; /* start with s as some value in a */
  for (i=0; i<m; i++) {
    for (j=0; j<n; j++) {
      if (a[i][j]<s) s=a[i][j];
    }
  }
  s-=1.; /* this makes sure the value ends up positive */
  for (i=0; i<m; i++) {
    for (j=0; j<n; j++) {
      a[i][j]-=s;
    }
  }

/*
 * Finish setting up the simplex problem.
 */

  rin=(int *)malloc(sizeof(*rin)*n);
  rout=(int *)malloc(sizeof(*rout)*m);
  cin=(int *)malloc(sizeof(*cin)*m);
  cout=(int *)malloc(sizeof(*cout)*n);
  for (i=0; i<m; i++) { 
    a[i][n]=1.;
    rout[i]=i;
    cin[i]= -1;
    }
  for (j=0; j<n; j++) {
    a[m][j]= -1.;
    cout[j]=j;
    rin[j]= -1;
  }
  *v=1.; /* Value of the game times the sum of the strategies */
  a[m][n]=0.;

/*
 * Simplex loop
 */

  while (1) { /* repeat until no suitable pivot is found */

/*
 * Get pivot
 */

    p= -1.; /* Mark that no pivot has been found yet */
    for (j=0; j<n; j++) { /* Take the largest of the column values */
      if (a[m][j]<0.) { /* Will a pivot in this column help? */
        cp= -2.; /* Mark that no pivot has been found in this column */
        for (i=0; i<m; i++) { /* Take the smallest value in each column */
          if (a[i][j]>1.e-6) { /* Suitable pivots must be positive */
            pc= -a[m][j]*a[i][n]/a[i][j];
            if (pc<cp||cp== -2.) {
              cp=pc;
              icp=i;
            }
          }
        }
        if (cp>p) {
          p=cp;
          ip=icp;
          jp=j;
        }
      }
    }

/*
 * If no pivot has been found the game is solved.
 */

    if (p== -1.) break;

/*
 * Perform the pivot operation.
 */

    for (i=0; i<=m; i++) { /* Be sure to include extra row */
      if (i!=ip) { /* Pivot row is unchanged except for the pivot */
	for (j=0; j<=n; j++) { /* Be sure to include extra column */
	  if (j!=jp) { /* Negate pivot column (not the pivot) later */
	    a[i][j]=(a[i][j]*a[ip][jp]-a[ip][j]*a[i][jp])/(*v);
	  }
	}
	a[i][jp]= -a[i][jp]; /* Now it's safe to change pivot column */
      }
    }
    p=a[ip][jp]; /* Don't update pivot until after we are done using it */
    a[ip][jp]= *v; /* Switch pivot value with the value of the game */
    *v=p;
    i=rin[jp]; /* Keep track of active strategies */
    rin[jp]=rout[ip];
    rout[ip]=i;
    j=cin[ip];
    cin[ip]=cout[jp];
    cout[jp]=j;
  }

/*
 * Copy over results.
 */

  *v=(*v)/a[m][n]+s;
  for (i=0; i<m; i++) r[i]=0; /* Set all row strategies to zero */
  for (j=0; j<n; j++) { /* Then copy over the values of the active strategies */
    if (rin[j]>=0) r[rin[j]]=a[m][j]/a[m][n];
  }
  for (j=0; j<n; j++) c[j]=0; /* Set all column strategies to zero */
  for (i=0; i<m; i++) { /* Then copy over the values of the active strategies */
    if (cin[i]>=0) c[cin[i]]=a[i][n]/a[m][n];
  }
  free(rin); /* don't let storage accumalate */
  free(rout);
  free(cin);
  free(cout);
}

/*
 * Function to solve AX=B by Guasian elimination.
 *
 * The rank of A is returned so that it is possible to tell when A is
 * singular and X is not returned.
 *
 * Parameters:
 * n      The size of the problem (input)
 * a      n x n matrix (input/scratch)
 * b      n vector (input/scratch)
 * x      n vector (output)
 */

int lineqs(n,a,b,x)

int n;
float **a;
float *b;
float *x;

{
int *isw; /* Table used to keep track of row switches */
int *jsw; /* Table used to keep track of column switches */
int jp,kp; /* Pivot location */
float pivot; /* Pivot value */
int i,j,k; /* Scratch */
float e; /* Scratch */

/*
 * Initialize row and column switch tables.
 */

  isw=(int *)malloc(sizeof(*isw)*n);
  jsw=(int *)malloc(sizeof(*jsw)*n);
  for (i=0; i<n; i++) {
    isw[i]=i;
    jsw[i]=i;
  }

/*
 * Triangularize a.
 */

  for (i=0; i<n; i++) {

/*
 * Use maximal element as pivot.
 */

    pivot=0.;
    for (j=i; j<n; j++) {
      for (k=i; k<n; k++) {
	e=a[isw[j]][jsw[k]];
	if (e<0.) e= -e;
	if (e>pivot) {
	  pivot=e;
	  jp=j;
	  kp=k;
	}
      }
    }

/*
 * If the pivot is essentially zero then a is singular. If so we can't
 * solve the problem and should just return the rank of a.
 */

    if (pivot<1e-5) {
      free(isw); /* don't let storage accumalate */
      free(jsw);
      return i;
    }

/*
 * Switch pivot row and column with the ith row and column.
 */

    j=isw[i];
    isw[i]=isw[jp];
    isw[jp]=j;
    k=jsw[i];
    jsw[i]=jsw[kp];
    jsw[kp]=k;

/*
 * Eliminate the remainder of the column.
 */

    pivot=a[isw[i]][jsw[i]];
    for (j=i+1; j<n; j++) {
      e= -a[isw[j]][jsw[i]]/pivot;
      for (k=i+1; k<n; k++) {
	a[isw[j]][jsw[k]]+=e*a[isw[i]][jsw[k]];
      }
      b[isw[j]]+=e*b[isw[i]];
    }
  }

/*
 * Back substitute to obtain x.
 */

  for (i=n-1; i>=0; i--) {
    x[jsw[i]]=b[isw[i]];
    for (j=n-1; j>i; j--) {
      x[jsw[i]]-=a[isw[i]][jsw[j]]*x[jsw[j]];
    }
    x[jsw[i]]/=a[isw[i]][jsw[i]];
  }
  free(isw); /* don't let storage accumalate */
  free(jsw);
  return n; /* matrix is not singular, so return the full rank */
}

This page is maintained by Bruno Wolff III on wolff.to.