/********************************************************************/
/*                                                                  */
/*  Module        : SolveEquations                                  */
/*                                                                  */
/*  Version       : 3.6                                             */
/*  Last revision : 03/28/94 20:14:24                               */
/*                                                                  */
/*  Description :                                                   */
/*     This module supplies a function for solving linear           */
/*     inhomogenous equations over F3. A special solution and a     */
/*     complete system of fundamental solutions are returned.       */
/*                                                                  */
/*  Functions supplied :                                            */
/*     int solve_equations ();                                      */
/*                                                                  */
/********************************************************************/

#include "aglobals.h"
#include "fdecla.h"
# include	"storage.h"

extern int prime;

char matrix[YMAX][XMAX];
int  x_dim, y_dim, yh_dim;
VEC  absolut, inhom;
VEC  fsolution[XMAX];

static char indicator[XMAX];
static long cxdim, cydim;

void showm ( void )
{
	long i, j;
      for ( i = 0; i < y_dim; i++ ) {
      	for ( j = 0; j < x_dim; j++ )
      		printf ( "%1d", matrix[i][j] );
      	printf ( "\n" ); 
      }
      puts ( "absolut" );
      write_vector ( absolut, y_dim );
}

/* internal functions */

/* F2 routines */

/*
void zero2_col ( row, col )
long row, col;
{
      register long i = y_dim;
	 register int wid;
	 register VEC r = matrix[row];
	 register VEC s;
	 
      while ( i-- ) {
            if ( i != row ) {
                  if ( matrix[i][col] != 0 ) {
			   		s = matrix[i];
			   		for ( wid = (int)col+1; wid--; )
						s[wid] ^= r[wid];
                       	absolut[i] ^= absolut[row];
                  }
            }
      }
} */

void zero2_col ( row, col )
long row, col;
{
 	register long i = y_dim;
	while ( i ) {
		i--;
		if ( i != row ) {
			if ( matrix[i][col] != 0 ) {
				add2_vector ( matrix[row], matrix[i], (int)col+1 );
				absolut[i] ^= absolut[row];
			}
		}
	}
}

void zeroh2_col ( row, col, fend )
long row, col;
int fend;
{
      register long i = fend;
      do {
            if ( i != row ) {
                  if ( matrix[i][col] != 0 )
                        add2_vector ( matrix[row], matrix[i], (int)cxdim );
            }
      } while ( ++i < cydim );
}

void zeroe2_col ( row, col, cstart )
long row, col;
int cstart;
{
	register long i = cydim;
	while ( i-- ) {
		if ( i != row ) {
			if ( matrix[i][col] != 0 )
				add2_vector ( matrix[row], matrix[i],
					(int)(cxdim+cstart) );
		}
	}
}

int gauss2_eliminate()
{
      register long ix;
      register long iy = 0;
      register char value;
      int rank = 0;
      int solvable = TRUE;
      
      zero_vector ( indicator, x_dim );
      while ( iy < y_dim && solvable ) {
            ix = x_dim - 1;
            while ( ( value = matrix[iy][ix] ) == 0 && ix ) ix--;
            if ( value ) {
                  rank++;
                  indicator[ix] = TRUE;
                  zero2_col ( iy, ix );
            }
            else
            	solvable = !absolut[iy];
		iy++;
      }
      if ( solvable )
	      return ( rank );
	else
		return ( -1 );
}

/* F3 routines */

void
zero3_col ( row, col )
long row, col;
{
      register long i = y_dim;
      register char x;
	while ( i--) {
            if ( i != row ) {
                  switch ( matrix[i][col] ) {
                        case 1:
                              x = absolut[i] + 3;
                              subb3_vector ( matrix[row], matrix[i], (int)col+1 );
                             	x -= absolut[row];
                             	absolut[i] = x > 2 ? x - 3 : x;
                              break;
                        case 2:
                              x = absolut[i];
                              add3_vector ( matrix[row], matrix[i], (int)col+1 );
                             	x += absolut[row];
                             	absolut[i] = x > 2 ? x - 3 : x;
                  }
            }
      }
}

void
zeroh3_col ( row, col, fend )
long row, col;
int fend;
{
      register long i = fend;
      do {
            if ( i != row ) {
                  switch ( matrix[i][col] ) {
                        case 1:
                              subb3_vector ( matrix[row], matrix[i], (int)cxdim );
                              break;
                        case 2:
                              add3_vector ( matrix[row], matrix[i], (int)cxdim );
                  }
            }
      } while ( ++i < cydim );
}

void zeroe3_col ( row, col, cstart )
long row, col;
int cstart;
{
	register long i = cydim;
	while ( i-- ) {
		if ( i != row ) {
			switch ( matrix[i][col] ) {
				case 1:
					subb3_vector ( matrix[row], matrix[i], (int)(cxdim+cstart) );
					break;
				case 2:
					add3_vector ( matrix[row], matrix[i], (int)(cxdim+cstart) );
			}
		}
	}
}

int gauss3_eliminate()
{
      register long ix;
      register long iy = 0;
      register char x, value;
      int rank = 0;
      int solvable = TRUE;
      
      zero_vector ( indicator, x_dim );
      while ( iy < y_dim && solvable ) {
            ix = x_dim - 1;
            while ( ( value = matrix[iy][ix] ) == 0 && ix ) ix--;
            switch ( value ) {
                  case 2:
                        smul3_vector ( 2, matrix[iy], x_dim );
                        rank++;
                        indicator[ix] = TRUE;
                        x = absolut[iy] << 1;
                        absolut[iy] = x > 3 ? 1 : x;
                        zero3_col ( iy, ix );
                        break;
                  case 1:
                        rank++;
                        indicator[ix] = TRUE;
                        zero3_col ( iy, ix );
                        break;
                  default:
                  	solvable = !absolut[iy];
            }
		iy++;
      }
      if ( solvable )
	      return ( rank );
	else
		return ( -1 );
}

/* Fp routines */

void
zerop_col ( row, col )
long row, col;
{
      register long i = y_dim;
      register char x, y;
	while ( i--) {
            if ( i != row ) {
                  if ( (x = matrix[i][col]) != 0 ) {
				sub_mult ( x, matrix[row], matrix[i], (int)col+1 );
	                  y = fp_mul ( x, absolut[row] );
	                  x = absolut[i] + prime - y;
                       	absolut[i] = x >= prime ? x - prime : x;
                  }
            }
      }
}

void
zerohp_col ( row, col, fend )
long row, col;
int fend;
{
      register long i = fend;
      register char x;
      do {
            if ( i != row )
                  if ( (x = matrix[i][col]) != 0 ) 
				sub_mult ( x, matrix[row], matrix[i], (int)col+1 );
      } while ( ++i < cydim );
}

void zeroep_col ( row, col, cstart )
long row, col;
int cstart;
{
	register long i = cydim;
	register char x;
	while ( i-- ) {
		if ( i != row ) {
			if ( (x = matrix[i][col]) != 0 ) 
				sub_mult ( x, matrix[row], matrix[i], (int)(cxdim+cstart) );
		}
	}
}

int gaussp_eliminate()
{
      register long ix;
      register long iy = 0;
      register char value;
      int rank = 0;
      int solvable = TRUE;
      
      zero_vector ( indicator, x_dim );
      while ( iy < y_dim && solvable ) {
            ix = x_dim - 1;
            while ( ( value = matrix[iy][ix] ) == 0 && ix ) ix--;
            if ( value != 0 ) {
            	value = fp_inv ( value );
	            smulp_vector ( value, matrix[iy], x_dim );
                  rank++;
                  indicator[ix] = TRUE;
                  absolut[iy] = fp_mul ( value, absolut[iy] );
                  zerop_col ( iy, ix );
		}
		else
                	solvable = !absolut[iy];
		iy++;
      }
      if ( solvable )
	      return ( rank );
	else
		return ( -1 );
}

/* common routines */

int fundamental_solutions()
{
      register long ix, iy;
      int solvable = TRUE;
      register char value;
      long j;

      for ( ix = 0; ix < x_dim; ix++ ) {
            if ( !indicator[ix] ) {
                  fsolution[ix] = ALLOCATE ( x_dim );
                  zero_vector ( fsolution[ix], x_dim );
                  *(fsolution[ix]+ix) = prime-1;
            }
            else
                  fsolution[ix] = NIL;
      }
      iy = 0;
      while ( iy < y_dim && solvable ) {
            ix = x_dim - 1;
            while ( ( value = matrix[iy][ix] ) == 0 && ix ) ix--;
            if ( value )
                  inhom[ix] = absolut[iy];
            else
                  solvable = !absolut[iy];
            for ( j = ix-1; j > -1; j-- ) {
                  if ( ( value = matrix[iy][j] ) != 0 )
                        *(fsolution[j]+ix) = value;
            }
            iy++;
      }
      return ( solvable );
}

/* end of internal functions */

int solve_equations ( int x, int y )
{
      int rank;
      
	 x_dim = x;
	 y_dim = y;
      zero_vector ( inhom, x_dim );
      rank = GAUSS_ELIMINATE();
      if ( ( rank != -1 ) && fundamental_solutions() )
		return ( rank );
      else
      	return ( -1 );
}

int gauss_p_eliminate ( x, y )
int x, y;
{
	register long ix, iy, i;
	register char value;
	int rank = 0;

	cxdim = x;
	cydim = y;
	
	for ( i = 0; i < y; i++ ) {
		zero_vector ( matrix[i]+(long)x, y );
		matrix[(long)i][(long)(x+i)] = 1;
	}

	for ( iy = 0; iy < y; iy++ ) {
		ix = x - 1;
		while ( ( value = matrix[iy][ix] ) == 0 && ix ) ix--;
		if ( value  ) {
			rank++;
			if ( value != 1 )
				SMUL_VECTOR ( fp_inv ( value ), matrix[iy], (int)(x+y) );
			ZEROE_COL ( iy, ix, y );
		}
	}
	return ( rank );
}

int complement ( c_start, cx_dim, cy_dim )
int c_start, cx_dim, cy_dim;
/*	Let {matrix[c_start],..,matrix[cy_dim-1]} be a generating set for a
	vector space V. Let S := {matrix[0],...,matrix[c_start-1]} be a subset
	of V. Then complement computes a basis for V/<S> (fsolution) and
	returns the dimension of this space (c_dim).
*/
{
      register long ix, iy, i;
      register char value;
      int c_dim = 0;

      cxdim = cx_dim;
      cydim = cy_dim;
      for ( iy = 0; iy < cy_dim; iy++ ) {
            ix = cxdim - 1;
            while ( ( value = matrix[iy][ix] ) == 0 && ix ) ix--;
            if ( value  ) {
            	if ( value != 1 )
            		SMUL_VECTOR ( fp_inv ( value ), matrix[iy], (int)cxdim );
                  if ( iy < c_start )
                        ZEROH_COL ( iy, ix, 0 );
                  else
                        ZEROH_COL ( iy, ix, c_start );
            }
      }
      for ( i = c_start; i < cy_dim; i++ ) {
            ix = cx_dim - 1;
            while ( ( value = matrix[i][ix] ) == 0 && ix ) ix--;
            if ( value ) {
            	fsolution[c_dim] = ALLOCATE ( cx_dim );
                  copy_vector ( matrix[i], fsolution[c_dim++], cx_dim );
		}
      }
      return ( c_dim );
}

/* end of module solveequations */