/* an_eqn.c -- a specific equation to solve with solveqn.c
 */

#include <stdio.h>
#include <stdlib.h>


// Find x such that A * x + b = 0.  Modifies A and b.
void solve( double **A, double *x, double *b, int size );


/* input is, e.g.,
   3, { "a", "b", "c" }, { 1, 2, 1, 3,  4, 5, 1, 6,   7, 8, 1, 9 }
   to mean
         1a + 2b + c = 3
         4a + 5b + c = 6
         7a + 8b + c = 9
   The last line in each row becomes negated in the constant array.
*/
void load( int size, double *coefs, double *copy, 
           double **A, double *x, double *b ) {
  int i, j, k;

  for( i = 0;  i < size;  i++ ) {
    for( j = 0;  j < size;  j++ ) {
      k = i * ( size + 1 ) + j;
      copy[k] = coefs[k];
    }
    k = i * ( size + 1 );
    A[i] = &copy[ k ];
    b[i] = -coefs[ k + size ];
  }
}


void load_solve_print(  int size, char **names, double *coefs ) {
  double *copy;
  double **A, *x, *b;

  double t;
  int i, j;

  copy = (double *) malloc( sizeof(double) * size * (size + 1) );
  A = (double **) malloc( sizeof(double *) * size );
  x = (double *) malloc( sizeof(double) * size );
  b = (double *) malloc( sizeof(double) * size );
  load( size, coefs, copy, A, x, b );
  solve( A, x, b, size );
  load( size, coefs, copy, A, x, b );
  for( i = 0;  i < size;  i++ ) {
    if( x[i] >= 0 ) {
      printf( "+ %g %s ", +x[i], names[i] );
    } else {
      printf( "- %g %s ", -x[i], names[i] );
    }
  }
  printf( "\n\n" );

  for( i = 0;  i < size;  i++ ) {
    t = 0.0;
    for( j = 0;  j < size;  j++ ) {
      if( A[i][j] == 1.0 ) {
	if( x[j] < 0 ) {
	  printf( "- %g. ", -x[j] );
	} else {
	  printf( "+ %g. ", +x[j] );
	}
      } else {
	printf( "+ %g*%g. ", A[i][j], x[j] );
      }
      t += A[i][j] * x[j];
    }
    printf( "= %g\n", t );
  }
  free( (void *) copy );
  free( (void *) A );
  free( (void *) x );
  free( (void *) b );
}


char *prob1_names [] = { "i", "(i%6)", "" }; // "" == constant's name
double prob1_coefs [] = {
  8, 2, 1, 30 - 8,
  9, 3, 1, 60 - 9,
  14, 2, 1, 150 - 14
};


char *prob2_names [] = { "p", "i", "(i%6)", "(i%30)", "" }; // "" == constant's name
double prob2_coefs [] = {
  40, 38, 2, 8, 1, 1 * 30,
  41, 39, 3, 9, 1, 2 * 30,
  46, 44, 2, 14, 1, 5 * 30,
  70, 68, 2, 8, 1, 17 * 30,
  512 + 40, 38, 2, 8, 1, 257 * 30
};



main() {
  load_solve_print( 3, prob1_names, prob1_coefs );
  printf( "\n\n" );
  load_solve_print( 5, prob2_names, prob2_coefs );
  printf( "\n\n" );
}