-- LinearSolver.mesa, a definition of a Cedar module for solving linear equations and inequations by gaussian elimination and the simplex method.  Part of the Juno 2 graphics programming system.  Coded by Greg Nelson, November 1983.
-- Last Edited by: Gnelson, December 11, 1983 10:45 pm

LinearSolver: DEFINITIONS = BEGIN

UnknownRec: TYPE = RECORD [r, name: REF, i, j: INT, b: BOOL]; 
Unknown: TYPE = REF UnknownRec; 
NewUnknown: PROC [name: REF] RETURNS [Unknown]; 

-- The client of LinearSolver creates (with NewUnknown) objects of type LinearSolver.Unknown, which are references to records.  A LinearSolver.Unknown represents an unknown to be solved for.  Clients should ignore the components of 

Push, Pop: PROC; 
Assert: PROC[c: Constraint]; 
Satisfiable: PROC RETURNS [BOOL]; 
Solution: PROC[v: Unknown] -- to be called only if Satisfiable[] -- RETURNS [REAL]; 

-- The module represents a stack of constraints, initially empty.  The client pushes a constraint onto the stack with Assert.  The effect of Push is to mark the current depth of the stack; the effect of Pop is to remove constraints from the stack until the last marked depth is restored.  Satisfiable[] returns true if the current constraints can be satisfied.  In this case, a solution to the constraints can be obtained by calling Solution[v] for each Unknown v.  Note that this solution may not be unique.  If you are interested in getting some particular solution instead of any solution, then you have to read about pivoting, below.

Constraint: TYPE = LIST OF REF -- to LinearMonomial --; 
LinearMonomial: TYPE = RECORD [coefficient: REAL, unknown: Unknown]; 
unity: Unknown;

-- "unity" is an Unknown representing one unit.  It exists before any Unknowns are created.  Solution[unity] = 1 invariantly.  The constraint represented by a list of linear monomials is that the some of them is zero.  Thus the constraint x+x=4 would be asserted Assert[[[2, x], [-4, unity]]], which asserts that 2*x+(-4)*unity = 0, or that x+x=4.

-- Example.  A madmen obsessed by a fear that 2 + 2 = 4 will fail to be true in the near future decides to write a program to keep checking and alert him as soon as it fails.  Since efficiency is not crucial, he chooses not to compute 2 + 2 and compare it with 4, but to solve x + x = 4 and compare the solution value of x with 2.  He could do so with the following code.

-- Madman's code.  
-- BEGIN OPEN LinearSolver; 
--   x: Unknown;  
--   WHILE TRUE DO 
--     Push[]; 
--     x := NewUnknown[];
--     Assert[LIST[[2, x], [-4, unity]]]; 
--     IF NOT Satisfiable[] OR Solution[x] # 2 THEN SIGNAL Madman; 
--     Pop[] 
--   ENDLOOP 
-- END.

Init: PROC;  -- Client must call it once before using the solver.

END. -- of LinearSolver.mesa