/*
 * 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 */
}