[alt.sources] lin: solve systems of linear equations

joe@dayton.UUCP (Joseph P. Larson) (01/05/90)

Okay.  Included is "lin", my rather small, cute program that solves
systems of n linear equations in m variables.  This version has been
changed since my very recent posting in that it solves systems with more
than one solution.  Of course, that means it'll say "all values" for
certain variables.  For instance, the solution for:

	1*X0      + 3*X2 = 4
		 2*X1 + 1*X2 = -2

will be:

	X0:  4 - 3*X2
	X1:  -1 - 0.5*X2
	X2:  all values

This program may be handy for people who are about to study Linear
Algebra.  However, I would not recommend you use lin as a substitute for
doing your problem sets by hand.   Of course, you *can* use it to check
your answers...

To create an executable, save to lin.c, edit off the header and the .sig
at the end, and type "make lin".

The program isn't very nice about how it displays floating point values.
Unfortunately, printf() is pretty stupid when you use the "%f" option and
I didn't feel like writing my own.  Thus, I just use "%8.4f".  If this
isn't enough precision for you -- change it or write your own routine to
convert floating point values to strings.

If you find this program helpful, drop me a line and pat me on the back.
It's nice to know sometimes that you aren't feeding an empty software pit.

-Joe

/*
 * lin: solve systems of linear equations.  This program is public domain.
 *
 * Description:
 *
 *		lin solves systems of linear equations.
 *
 *		See the code or the "usage" line for any help you need running it.
 *
 * Author:
 *
 *		Joe Larson
 *		(currently joe@dayton.dhdsc.mn.org)
 *
 * Revision History
 *
 *		12-89	jpl.	Created.
 *		 1-90	jpl.	Handle systems w/ infinite solutions.
 */

#include <stdio.h>
#include <ctype.h>

#define		MAXVAR	20
#define		MAXEQN	20

double strtod();
short	vcount, ecount, solved[MAXVAR];
double	matrix[MAXEQN][MAXVAR+1];
char	showstep = 0;

main(argc,argv)
int argc;
char **argv;
{
	char	*c,
			line[132],
			x;
	short	i,j, evmin;

	for (i = 1; i < argc; i++)
	{
		if (!strcmp(argv[i], "-s")) showstep = 1;
		else usage(argv[0]);
	}
	printf("This program solves systems of linear equations involving\n");
	printf("n variables and m unknowns.\n");
	while(1)
	{
		while (1)
		{
			printf("Number of variables: "); fflush(stdout);
			gets(line);
			vcount = atoi(line);
			if (vcount > MAXVAR)
			{
				printf("Max variables = %d.  ", MAXVAR);
				continue;
			}
			if (vcount) break;
			exit(0);
		}
		while (1)
		{
			printf("Number of equations: "); fflush(stdout);
			gets(line);
			ecount = atoi(line);
			if (ecount > MAXEQN)
			{
				printf("Max equations = %d.  ", MAXEQN);
				continue;
			}
			if (ecount) break;
			exit(0);
		}
		evmin = (ecount < vcount) ? ecount : vcount;

		printf("Equations are of the form: \n");
		for (i = 0; i < vcount; i++)
		{
			if (i) printf("+ ");
			printf(" c%d * x%d ", i, i);
		}
		printf("= result\n");
		printf("For each equation, enter the %d constants plus the result,\n",
				vcount);
		printf("separated by spaces.\n");

		for (i = 0; i < ecount; i++)
		{
			printf("\nEquation %d\n", i+1);
			gets(line);
			for (c = line, j = 0; *c && j < vcount; j++)
				matrix[i][j] = strtod(c, &c);
			matrix[i][vcount] = strtod(c, &c);
		}

		printf("The system: \n");
		show();

	/*
	 * Now solve the equation in the following manner:
	 *
	 *		for each variable:
	 *				find an equation with a 1 for constant.
	 *					(if non, find one with a non-zero for constant
	 *					and multiply to make it a 1.)
	 *				add multiples of this equiation to others as needed
	 *					to zero-out their constants.
	 */

		for (i = 0; i < evmin; i++)
		{
			/* Find an equation to place here */
			x = 0;
			for (j = i; j < ecount; j++)
				if (matrix[j][i] == 1.0)
				{
					swap(i, j);
					x = 1;
					break;
				}
			if (!x)
				for (j = i; j < ecount; j++)
					if (matrix[j][i] != 0.0)
					{
						swap(i,j);
						x = 1;
						mult(i, 1.0 / matrix[i][i]);
						break;
					}
			if (!x)
			{
				/* This variable doesn't matter. */
				continue;
			}

			/*
			 * Now -- zero out everyone else.
			 */
			for (j = 0; j < ecount; j++)
			{
				if (j == i) continue;
				if (matrix[j][i] == 0.0) continue;
				add(j, i, -matrix[j][i]);
			}
		}

		/*
		 * Now -- solve for each of the variables.
		 */
		solvecomplex();

		printf("The system: \n");
		show();
	}
}

solvecomplex()
/*
 *
 */
{
	short i, j, k, lastj;

	for (i = 0; i < ecount; )
	{
		for (j = 0; j < vcount; j++)
			if (matrix[i][j]) break;
		if (j == vcount)
		{
			if (matrix[i][j])	/* No solutions. */
			{
				printf("The set has no real-number solutions.\n");
				return;
			}
			swap(i,ecount-1);
			ecount--;
			continue;
		}
		i++;
	}
	if (!ecount)
	{
		printf("System works for all values.\n");
		return;
	}
	lastj = -1;
	for (i = 0; i < ecount; i++)
	{
		for (j = 0; j < vcount; j++)
			if (matrix[i][j] != 0.0) break;
		if (j != i)
			for (lastj++; lastj < j; lastj++)
				printf("X(%d): all values.\n");
		lastj = j;
		printf("X(%d): %8.4f", j, matrix[i][vcount]);
		for (k = j+1; k < vcount; k++)
			if (matrix[i][k] != 0.0)
				printf(" - [%8.4f * X(%d)]", matrix[i][k], k);
		printf("\n");
	}
	for (i = lastj + 1; i < vcount; i++)
		printf("X(%d): all values\n", i);
}

swap(i,j)
short i,j;
/*
 * Swap two rows of the matrix.
 */
{
	double	temp;
	short	x;

	if (i == j) return;
	for (x = 0; x <= vcount; x++)
	{
		temp = matrix[i][x];
		matrix[i][x] = matrix[j][x];
		matrix[j][x] = temp;
	}
	if (showstep)
	{
		printf("Swap row %d and %d leaving:\n", i, j);
		show();
	}
}

mult(i, f)
short i;
double f;
/*
 * Multiply row i by constant f.
 */
{
	short x;

	for (x = 0; x <= vcount; x++)
		matrix[i][x] *= f;
	if (showstep)
	{
		printf("Multiply row %d by %f leaving:\n", i, f);
		show();
	}
}

add(i,j,f)
short i,j;
double f;
/*
 * Add to row i row (j*f).
 */
{
	short x;

	for (x = 0; x <= vcount; x++)
		matrix[i][x] += matrix[j][x] * f;
	if (showstep)
	{
		printf("Add (%f times row %d) to row %d:\n", f, j, i);
		show();
	}
}

show()
/*
 * Show the equation
 */
{
	short i,j;

	for (i = 0; i < ecount; i++)
	{
		for (j = 0; j < vcount; j++)
			printf(" %8.4f", matrix[i][j]);
		printf("   %8.4f\n", matrix[i][j]);
	}
}

usage(arg)
char *arg;
{
	fprintf(stderr, "\
Usage: %d [-s]\n\
where -s selects display of steps used to solve the set of equations.\n", arg);
	
	exit(0);
}
-- 
UUCP: rutgers!dayton!joe   (Picts 1-16 are   DHDSC - Joe Larson/MIS 1060
ATT : (612) 375-3537       now ready.)       700 on the Mall, Mpls, Mn. 55402
Newsgroups: alt.source
Subject: lin: solve systems of linear equations (new version)
Reply-To: joe@dayton.UUCP (Joseph P. Larson)
Distribution: alt
Organization: Dayton-Hudson Dept. Store Co.


Okay.  Included is "lin", my rather small, cute program that solves
systems of n linear equations in m variables.  This version has been
changed since a very recent posting in that it solves systems with more
than one solution.  Of course, that means it'll say "all values" for
certain variables.  For instance, the solution for:

	1*X0      + 3*X2 = 4
		 2*X1 + 1*X2 = -2

will be:

	X0:  4 - 3*X2
	X1:  -1 - 0.5*X2
	X2:  all values

This program may be handy for people who are about to study Linear
Algebra.  However, I would not recommend you use lin as a substitute for
doing your problem sets by hand.   Of course, you *can* use it to check
your answers...

To create an executable, save to lin.c, edit off the header and the .sig
at the end, and type "make lin".

The program isn't very nice about how it displays floating point values.
Unfortunately, printf() is pretty stupid when you use the "%f" option and
I didn't feel like writing my own.  Thus, I just use "%8.4f".  If this
isn't enough precision for you -- change it or write your own routine to
convert floating point values to strings.

If you find this program helpful, drop me a line and pat my on the back.
It's nice to know sometimes that you aren't feeding an empty software pit.

-Joe

/*
 * lin: solve systems of linear equations.  This program is public domain.
 *
 * Description:
 *
 *		lin solves systems of linear equations.
 *
 *		See the code or the "usage" line for any help you need running it.
 *
 * Author:
 *
 *		Joe Larson
 *		(currently joe@dayton.dhdsc.mn.org)
 *
 * Revision History
 *
 *		12-89	jpl.	Created.
 *		 1-90	jpl.	Handle systems w/ infinite solutions.
 */

#include <stdio.h>
#include <ctype.h>

#define		MAXVAR	20
#define		MAXEQN	20

double strtod();
short	vcount, ecount, solved[MAXVAR];
double	matrix[MAXEQN][MAXVAR+1];
char	showstep = 0;

main(argc,argv)
int argc;
char **argv;
{
	char	*c,
			line[132],
			x;
	short	i,j, evmin;

	for (i = 1; i < argc; i++)
	{
		if (!strcmp(argv[i], "-s")) showstep = 1;
		else usage(argv[0]);
	}
	printf("This program solves systems of linear equations involving\n");
	printf("n variables and m unknowns.\n");
	while(1)
	{
		while (1)
		{
			printf("Number of variables: "); fflush(stdout);
			gets(line);
			vcount = atoi(line);
			if (vcount > MAXVAR)
			{
				printf("Max variables = %d.  ", MAXVAR);
				continue;
			}
			if (vcount) break;
			exit(0);
		}
		while (1)
		{
			printf("Number of equations: "); fflush(stdout);
			gets(line);
			ecount = atoi(line);
			if (ecount > MAXEQN)
			{
				printf("Max equations = %d.  ", MAXEQN);
				continue;
			}
			if (ecount) break;
			exit(0);
		}
		evmin = (ecount < vcount) ? ecount : vcount;

		printf("Equations are of the form: \n");
		for (i = 0; i < vcount; i++)
		{
			if (i) printf("+ ");
			printf(" c%d * x%d ", i, i);
		}
		printf("= result\n");
		printf("For each equation, enter the %d constants plus the result,\n",
				vcount);
		printf("separated by spaces.\n");

		for (i = 0; i < ecount; i++)
		{
			printf("\nEquation %d\n", i+1);
			gets(line);
			for (c = line, j = 0; *c && j < vcount; j++)
				matrix[i][j] = strtod(c, &c);
			matrix[i][vcount] = strtod(c, &c);
		}

		printf("The system: \n");
		show();

	/*
	 * Now solve the equation in the following manner:
	 *
	 *		for each variable:
	 *				find an equation with a 1 for constant.
	 *					(if non, find one with a non-zero for constant
	 *					and multiply to make it a 1.)
	 *				add multiples of this equiation to others as needed
	 *					to zero-out their constants.
	 */

		for (i = 0; i < evmin; i++)
		{
			/* Find an equation to place here */
			x = 0;
			for (j = i; j < ecount; j++)
				if (matrix[j][i] == 1.0)
				{
					swap(i, j);
					x = 1;
					break;
				}
			if (!x)
				for (j = i; j < ecount; j++)
					if (matrix[j][i] != 0.0)
					{
						swap(i,j);
						x = 1;
						mult(i, 1.0 / matrix[i][i]);
						break;
					}
			if (!x)
			{
				/* This variable doesn't matter. */
				continue;
			}

			/*
			 * Now -- zero out everyone else.
			 */
			for (j = 0; j < ecount; j++)
			{
				if (j == i) continue;
				if (matrix[j][i] == 0.0) continue;
				add(j, i, -matrix[j][i]);
			}
		}

		/*
		 * Now -- solve for each of the variables.
		 */
		solvecomplex();

		printf("The system: \n");
		show();
	}
}

solvecomplex()
/*
 *
 */
{
	short i, j, k, lastj;

	for (i = 0; i < ecount; )
	{
		for (j = 0; j < vcount; j++)
			if (matrix[i][j]) break;
		if (j == vcount)
		{
			if (matrix[i][j])	/* No solutions. */
			{
				printf("The set has no real-number solutions.\n");
				return;
			}
			swap(i,ecount-1);
			ecount--;
			continue;
		}
		i++;
	}
	if (!ecount)
	{
		printf("System works for all values.\n");
		return;
	}
	lastj = -1;
	for (i = 0; i < ecount; i++)
	{
		for (j = 0; j < vcount; j++)
			if (matrix[i][j] != 0.0) break;
		if (j != i)
			for (lastj++; lastj < j; lastj++)
				printf("X(%d): all values.\n");
		lastj = j;
		printf("X(%d): %8.4f", j, matrix[i][vcount]);
		for (k = j+1; k < vcount; k++)
			if (matrix[i][k] != 0.0)
				printf(" - [%8.4f * X(%d)]", matrix[i][k], k);
		printf("\n");
	}
	for (i = lastj + 1; i < vcount; i++)
		printf("X(%d): all values\n", i);
}

swap(i,j)
short i,j;
/*
 * Swap two rows of the matrix.
 */
{
	double	temp;
	short	x;

	if (i == j) return;
	for (x = 0; x <= vcount; x++)
	{
		temp = matrix[i][x];
		matrix[i][x] = matrix[j][x];
		matrix[j][x] = temp;
	}
	if (showstep)
	{
		printf("Swap row %d and %d leaving:\n", i, j);
		show();
	}
}

mult(i, f)
short i;
double f;
/*
 * Multiply row i by constant f.
 */
{
	short x;

	for (x = 0; x <= vcount; x++)
		matrix[i][x] *= f;
	if (showstep)
	{
		printf("Multiply row %d by %f leaving:\n", i, f);
		show();
	}
}

add(i,j,f)
short i,j;
double f;
/*
 * Add to row i row (j*f).
 */
{
	short x;

	for (x = 0; x <= vcount; x++)
		matrix[i][x] += matrix[j][x] * f;
	if (showstep)
	{
		printf("Add (%f times row %d) to row %d:\n", f, j, i);
		show();
	}
}

show()
/*
 * Show the equation
 */
{
	short i,j;

	for (i = 0; i < ecount; i++)
	{
		for (j = 0; j < vcount; j++)
			printf(" %8.4f", matrix[i][j]);
		printf("   %8.4f\n", matrix[i][j]);
	}
}

usage(arg)
char *arg;
{
	fprintf(stderr, "\
Usage: %d [-s]\n\
where -s selects display of steps used to solve the set of equations.\n", arg);
	
	exit(0);
}
-- 
UUCP: rutgers!dayton!joe   (Picts 1-16 are   DHDSC - Joe Larson/MIS 1060
ATT : (612) 375-3537       now ready.)       700 on the Mall, Mpls, Mn. 55402