PROGRAM simequ;

{
A program to demonstrate simultaneous equation solving using the Gaussian elimination method, as
described in many Engineering Mathematics texts.

Works with as many unknowns as required, the maximum is set by the available memory.


IMPORTANT!!!!!
########################################
# Compile with Range checking disabled #
########################################



Initial coding				06.01.95 to 09.01.95
Modified for release		15.01.95

	Release delayed by a damn near untraceable bug

Fixed the bug!				24.02.95
Tweaking and twiddling	18.03.95
Ditto							22.10.95 (It's FINALLY ready for release :-)
}

{
-------------------------------------------------------------------------------------------------------
--- TYPES, CONSTANTS AND VARIABLES --- TYPES, CONSTANTS AND VARIABLES --- TYPES, CONSTANTS AND VARIABLE
-------------------------------------------------------------------------------------------------------
}

{
Uses the old trick for defining dynamic arrays of setting the range to zero, then pointing a
dynamically allocated block of memory to the array variable, and accessing the memory via
normal array accesses.  Hence the need to compile with range checking disabled.
}
	
TYPE
	coeffictype = extended;
	coefficsarray = ARRAY[0..0] of coeffictype;
	resultsarray = ARRAY[0..0] of coeffictype;
	
VAR
	coeffics : ^coefficsarray;
	results : ^resultsarray;
	unknowns : longint;
	
{
-------------------------------------------------------------------------------------------------------
--- PROCEDURES AND FUNCTIONS --- PROCEDURES AND FUNCTIONS --- PROCEDURES AND FUNCTIONS --- PROCEDURES A
-------------------------------------------------------------------------------------------------------
}

FUNCTION offset(x, y : longint) : longint;

{
Converts between 2-D matrix coordinates and 1-D array offset.  Required since the dynamic array
technique doesn't seem to deal with arrays other than 1-D ones, but it's easier to work with 2-D
arrays for Gaussian Elimination.  Therefore we need this layer of translation.

Incidentally, the "damn near untraceable bug" was in this function.  As you can imagine, a bug
here will cause all memory access to be incorrect, leading to random program failure depending
on the state of the rest of the system memory.
}

BEGIN
	IF y > 0 THEN
	BEGIN
		offset := ((x - 1) + ((y - 1) * (unknowns + 1)))
	END
	ELSE
	BEGIN
		offset := (x - 1)
	END
END; {offset}


PROCEDURE show_matrix;

{
Prints the matrix, with a vertical bar separating the matrix A from the column vector b.

By default this is not used, but since it may be of use, the code is here if you need it.
}

VAR
	loop_a : integer;
	loop_b : integer;
	
BEGIN
	writeln;
	FOR loop_a := 1 TO unknowns DO
	BEGIN
		FOR loop_b := 1 TO unknowns DO
		BEGIN
			write(coeffics^[offset(loop_b,loop_a)]:0:3,'  ')
		END;
		writeln(' | ',coeffics^[offset(unknowns+1,loop_a)]:0:3)
	END
END; {show_matrix}


PROCEDURE show_results;

{
Prints the results for all the (previously) unknown terms
}

VAR
	loop : integer;
	result : coeffictype;
	
BEGIN
	writeln;
	writeln('Results :');
	writeln;
	FOR loop := 1 TO unknowns DO
	BEGIN
		result := results^[offset(loop,0)];
		writeln('x',loop,'=',result:0:3)
	END
END; {show_results}

		
PROCEDURE enter_equations;

{
Prompts user for the number of unknowns to solve, and then for all the coefficients and equation
answers.
}

VAR
	loop_a : integer;
	loop_b : integer;
	coeffic : coeffictype;
	
BEGIN
	FOR loop_a := 1 TO unknowns DO
	BEGIN
		writeln('Equation ',loop_a);
		writeln;
		FOR loop_b := 1 TO unknowns DO
		BEGIN
			write('x',loop_b,' = ');
			readln(coeffic);
			coeffics^[offset(loop_b,loop_a)] := coeffic
		END;
		write('Result = ');
		readln(coeffic);
		coeffics^[offset(unknowns + 1,loop_a)] := coeffic;
		writeln
	END
END; {enter_equations}


PROCEDURE prettyprint;

{
Prints equations in a somewhat neatened up form, with unknowns represented as x1, x2 etc. etc.
}

VAR
	loop_a : integer;
	loop_b : integer;
	coeffic : coeffictype;
	sign : boolean;
	
BEGIN
	writeln;
	FOR loop_a := 1 TO unknowns DO
	BEGIN
		FOR loop_b := 1 TO unknowns DO
		BEGIN
			coeffic := coeffics^[offset(loop_b,loop_a)];
			sign := (coeffic - ABS(coeffic) = 0);
			IF (loop_b > 1) THEN
			BEGIN
				CASE sign OF
					true : write(' +');
					false : write(' -')
				END;
				write(ABS(coeffic):0:3,'x',loop_b)
			END
			ELSE
			BEGIN
				write(coeffic:0:3,'x',loop_b)
			END
		END;
		coeffic := coeffics^[offset(unknowns + 1,loop_a)];
		writeln('=',coeffic:0:3)
	END
END; {prettyprint}


PROCEDURE gaussian_elimination;

{
Performs the Gaussian Elimination method of obtaining a triangular matrix.  I'm not going to
describe it here, you'll have to look in a maths text book.  Any half decent one, especially
Engineering Maths texts, will describe it in detail.
}

VAR
	loop_a : integer;
	loop_b : integer;
	loop_c : integer;
	row : integer;
	tosub : coeffictype;
	numer : coeffictype;
	denom : coeffictype;
	
BEGIN
	row := 1;
	FOR loop_a := 2 TO unknowns DO
	BEGIN
		denom := coeffics^[offset(row,row)];
		FOR loop_b := row + 1 TO unknowns DO
		BEGIN
			numer := coeffics^[offset(row,loop_b)];
			FOR loop_c := row TO unknowns DO
			BEGIN
				tosub := coeffics^[offset(loop_c,loop_b)];
				tosub := tosub - (coeffics^[offset(loop_c,row)] * (numer / denom));
				coeffics^[offset(loop_c,loop_b)] := tosub
			END;
			tosub := coeffics^[offset(unknowns + 1,loop_b)];
			tosub := tosub - (coeffics^[offset(unknowns + 1,row)] * (numer / denom));
			coeffics^[offset(unknowns + 1,loop_b)] := tosub			
		END;
		inc(row)
	END
END; {gaussian_elimination}


PROCEDURE results_time;

{
Determines the numerical values for each of the unknowns, working backwards from the last
unknown, as required by the use of an upper triangular matrix.
}

VAR
	loop_a : integer;
	loop_b : integer;
	row : integer;
	subresult : coeffictype;
	result : coeffictype;

	res : coeffictype;
	mul : coeffictype;
	dvs : coeffictype;

BEGIN
	row := unknowns - 1;
	result := coeffics^[offset(unknowns + 1,unknowns)] / coeffics^[offset(unknowns,unknowns)];
	results^[offset(unknowns,0)] := result;
	FOR loop_a := 1 TO unknowns - 1 DO
	BEGIN
		subresult := coeffics^[offset(unknowns + 1,row)];
		FOR loop_b := row + 1 TO unknowns DO
		BEGIN
			subresult := subresult - (coeffics^[offset(loop_b,row)] * results^[offset(loop_b,0)])
		END;
		result := subresult / coeffics^[offset(row,row)];
		results^[offset(row,0)] := result;
		dec(row)
	END
END; {results_time}


FUNCTION reserve_memory : boolean;

{
Attempts to reserve enough heap memory to store the matrices. Checks the largest block on the
heap to see if it's large enough. For the maximum of 10000 unknowns, you'd require somewhere in 
the region of.......954MBytes. So, either you have a computer the size of Mount Everest to hold
all the memory chips, or you have a damn big virtual memory swap file. Remember, the memory
requirement of this program is proportional the the SQUARE of the number of unknowns.
}

VAR
	block : longint;
	
BEGIN
	block := (((unknowns + 1) * unknowns) + 1) * sizeof(coeffictype);
	IF (block < maxavail) THEN
	BEGIN
		block := ((unknowns + 1) * unknowns) * sizeof(coeffictype);
		getmem(coeffics,block);
		block := unknowns * sizeof(coeffictype);
		getmem(results,block);
		reserve_memory := true
	END
	ELSE
	BEGIN
		reserve_memory := false
	END
END; {reserve_memory}


PROCEDURE free_memory;

{
Being a good little programmer, I'll free the memory I used.....or at least I'll
ask the system to free it for me.  Whether the system actually frees it is another
matter, but that's no concern of mine.
}

VAR
	blocksize : longint;
	
BEGIN
	blocksize := ((unknowns + 1) * unknowns) * sizeof(coeffictype);
	freemem(coeffics,blocksize);
	blocksize := unknowns * sizeof(coeffictype);
	freemem(results,blocksize)
END; {free_memory)
				
{
-------------------------------------------------------------------------------------------------------
--- PROGRAM CORE --- PROGRAM CORE --- PROGRAM CORE --- PROGRAM CORE --- PROGRAM CORE --- PROGRAM CORE -
-------------------------------------------------------------------------------------------------------
}
	
BEGIN
	writeln;
	writeln('Simultaneous Equations Solver');
	writeln('Version 1.0 - Compilation Date 22.10.95');
	writeln;
	write('How many unknowns :');
	readln(unknowns);
	IF (unknowns > 1) AND (unknowns <= 10000) THEN
	BEGIN
		IF reserve_memory THEN
		BEGIN
			enter_equations;
			prettyprint;
			gaussian_elimination;
			{show_matrix;}
			results_time;
			show_results;
			free_memory
		END
		ELSE
		BEGIN
			writeln('Sorry, insufficient free memory.')
		END
	END
	ELSE
	BEGIN
		writeln('Invalid number of unknowns.')
	END;
	writeln
END.