DIRECTORY

  RealFns,
  HerculesStorage,
  HerculesAlgebra USING [Frame, Matrix, PointPtr, GetFrameMatrix, TransformPoint],
  HerculesSolver,
  HerculesGraphics,
  Real,
  IO;

HerculesSolverImpl: PROGRAM

  IMPORTS
     RealFns, HerculesAlgebra, HerculesGraphics, Real
  EXPORTS
     HerculesSolver  =

BEGIN

OPEN
  Alg: HerculesAlgebra,
  Stor: HerculesStorage,
  HerculesSolver,
  Gr: HerculesGraphics;

Frame: TYPE = Alg.Frame;

PointPtr: TYPE = Alg.PointPtr;

displayingMotions: PUBLIC BOOL _ FALSE;

loggingPivots: BOOL _ FALSE;

rowmax: INTEGER = 100;

colmax: INTEGER = 100;

Column: TYPE = ARRAY[1..rowmax] OF REAL;

Tableau: TYPE = ARRAY[1..rowmax] OF ARRAY[1..colmax] OF REAL;

a: REF Tableau _ NEW[Tableau];

b: REF Column _ NEW[Column];

Tag: TYPE = {x, y, eq, geq};

TableauLabel: TYPE = RECORD [point: PointPtr, tag: Tag];

rowhead: ARRAY[1..rowmax] OF TableauLabel;

colhead: ARRAY[1..colmax] OF TableauLabel;

TooManyConstraints: SIGNAL = CODE;

TooManyConstrainedPoints: SIGNAL = CODE;

Square: PROCEDURE[x:REAL] RETURNS [REAL] =  BEGIN RETURN [x*x] END;

Sqrt: PROCEDURE[x:REAL] RETURNS [REAL] =  BEGIN RETURN [RealFns.SqRt[x]] END;

mFrame: REF Alg.Matrix _ NEW [Alg.Matrix];

Instr: TYPE = RECORD -- one step of compiled constraints
  [op: OpCode, -- step type
   larg: Row, ln: NRows, -- rows containing first arg of step
   rarg: Row, rn: NRows, -- rows containing second arg of step, if any
   link: REF Instr];

prog: RECORD -- compiled constraints
  [first, last: REF Instr, -- first and last instructions. Nodes from last.link on are free.
   n: NRows, -- Number of rows of a in use when prog ends.
   m: NCols -- Number of columns of a used (i.e., variables entering into the constraints)
   ]

CompileExpr: PROC [e: Se, alist: Alist] RETURNS [k: NRows] =
  -- Compiles on expression, appending its code to prog.
  --    The compiled code should leave the numerical components of the
  --    expression, and their gradients, on the last k rows of a.

BEGIN
  ???
END;

CompileConstr: PROC [ce: Se, alist: Alist] RETURNS [k: NRows] =
  -- Compiles a constraint expression, appending its code to prog.
  --    The compiled code should leave the error functions of the
  --    constraint (with appropriate tags) and their gradients on the last k rows of a.

BEGIN
  ???
END;

Solve: PUBLIC PROC[points: LIST OF NumPtr, constrs: LIST OF Se]
  RETURNS [success, impossible: BOOL] =

BEGIN

MarkMobilePoints: PROC [bool: BOOL] =
  -- Mark/unmark given points as mobile (unless they are fixed)

BEGIN
  p: PointPtr _ points;
  UNTIL p = NIL DO
    p.mobile  _ bool AND NOT p.fixed;
    p _ p.link
   ENDLOOP
END;

DisplayMotions: PROC [points: PointList] =
  -- Draws on the screen the motions of solved points

BEGIN
  p: PointPtr _ points;
  UNTIL p = NIL DO
     IF p.oldx # Real.RoundLI[p.x] OR p.oldy # Real.RoundLI[p.y]
        THEN Gr.DrawEdge[p.oldx, p.oldy, p.x, p.y];
     p _ p.link
   ENDLOOP
END;

NewtonRaphson: PROC RETURNS [success, impossible, moved: BOOLEAN] =
   -- Does one Newton-Raphson step.
   -- Returns impossible=TRUE if it concludes it is impossible to satisfy the current
   --   constraints given the current positions of the fixed and tempfixed points.
   -- Returns success=TRUE if all constraints are satisfied.
   -- Returns moved = TRUE if the coordinates of one or more points,
   --   rounded to the nearest integer, have changed from their previous rounded values. 

BEGIN

n, m: INT _ 0; -- the tableau a has n relevant rows and m relevant columns.

EvalWithGrad: PROC[ex: Se] RETURNS [results: Results, k: INTEGER] =
  -- Evaluates a symbolic expression ex that returns a structure with
  --     k (zero or more) real-valued simple components (no ropes, funvals, atoms). 
  -- Then adds k new rows to a[][] and b[]. Stores in b[] the 
  --   components of the result, and in a[][] the partial derivatives of each component
  --   with respect ot all mobile variables entering into it.
BEGIN
   WITH ex SELECT FROM
      NIL => RETURN [0];
      ATOM => {evaluate the atom, enumerate NumCells n result, push each};
      NumCell => {must be constant; anyway, check it}
      RopeCell=> error
      FunPtr => error
      LIST OF Se =>
         evaluate op = first component
         with op select from
            +, - => evaluate each operand recursively, check if
                     same num. of components, add/subtract rows of a[][] and b[]
           *, / => evaluate each operand; one must be single-valued
                    proceed as before
           xcoord, ycoord => eval operand, select specified row
           leftparen => eval operand
           lefbrack => eval operand, make into list
           comma => eval both sides, concat
           etc
END;

EvalWithGrad: PROC[vex: Se] =
  -- vex is an expression that returns two real values (coordinates of some 2-vector) 
  -- This procedure adds two new rows to the matrix a. It evaluates vex, and
  --    saves the coordinates of the result in the two elements of b[].
  -- The partial derivatives of each coordinate of vex
  --    with respect to all mobile NumCells referenced by vex
  --    are stored in the corresponding rows of a[][].

BEGIN
   ???
END;

EvNum: PROC[nex: Se] =
  -- nex is an expression that returns one real value 
  -- This procedure evaluates nex, and adds one new row to the matrix a.
  -- The value of nex is stored in
  --    the b[] vector. The partial derivatives of nex
  --    with respect to all mobile NumCells referenced by nex
  --    are stored in the corresponding rows of a[][].

BEGIN
   ???
END;

EvRelVector: PROC[vex: Se, fex: Se] =
  -- vex is an expression that returns one two real values (coordinates of some point) 
  -- fex is an expression that returns a frame (one to three vectors) 
  -- This procedure adds two new rows to the matrix a.
  --  It computes u = (vex times the pseudo-inverse of the matrix of fex)
  --   and stores the coordinates of u in the two elements of b[].
  --  The partial derivatives of these coordinates 
  --    with respect to all mobile NumCells referenced by vex or fex
  --    are stored in the corresponding rows of a[][].

BEGIN
   EvVector[vex];
   ??? <evalute the frame, 
END;

EvConstraint: PROC[cex: Se, fex: Se] =
  -- cex is a constraint expression that returns one or two real values
  --   "error functions" for the constraint) 
  -- fex is an expression that returns a frame (one to three vectors), or NIL
  -- This procedure adds one or two new rows to the matrix a.
  -- The error functions are stored in one or two elements of b[].
  -- The partial derivatives of these error functions 
  --    with respect to all mobile NumCells referenced by cex or fex
  --    are stored in the corresponding rows of a[][].
  -- The tag of each row is set ot eh appropriate constraint type (eq, geq, leq).
BEGIN
   ???
END;

AddRow: PROC[i, j: [1..rowmax], coef: REAL] =
  -- adds to row i of the matrix the contents of row j multiplied by coef.
  -- The tag of each row is set ot eh appropriate constraint type (eq, geq, leq).
BEGIN
   b[i] _ b[i] + coef * b[j];
   FOR c: INTEGER [1..colmax] DO 
      a[i][c] _ a[i][c] + coef * a[j][c]
    ENDLOOP
END;

PopRows: PROC [k: INTEGER] =
  -- Deletes the last k rows of the matrix.
BEGIN
   n _ n-k
END;

Constraint: PROC[value: REAL, tag: Tag] =
  -- adds row to matrix for a new constraint whose current value is value.

BEGIN
   n _ n + 1; -- add new row
   IF n > rowmax THEN SIGNAL TooManyConstraints[];
   rowhead[n] _ [point: NIL, tag: tag];
   IF tag=eq AND ABS[value] > 0.01 OR tag=geq AND value <-0.01 THEN
      {success _ FALSE};
   FOR j: INT IN [1..m] DO a[n][j] _ 0.0 ENDLOOP;
   b[n] _ - value;
END;

Partial: PROC[point: PointPtr, xCoef, yCoef: REAL] =
  -- Adds to the matrix one or two new entries, corresponding to the partial derivatives
  --   of the most recently entered error function
  --   with respect to the x and y coordinates of the given point.

BEGIN
   IF NOT point.mobile THEN RETURN;
   IF ABS[xCoef] > .000001 THEN
     BEGIN

        IF point.xCol = 0 THEN -- add new column:
          {m _ m + 1;
           IF m > colmax THEN SIGNAL TooManyConstrainedPoints[];
           FOR i: INT IN [1..n] DO a[i][m] _ 0.0 ENDLOOP;
           point.xCol _ m;
           colhead[m] _ [point:point, tag: x]};
        a[n][point.xCol] _ a[n][point.xCol] + xCoef

     END;

   IF ABS[yCoef] > .000001 THEN
     BEGIN

        IF point.yCol = 0 THEN -- add new column:
          {m _ m + 1;
           IF m > colmax THEN SIGNAL TooManyConstrainedPoints[];
           FOR i: INT IN [1..n] DO a[i][m] _ 0.0 ENDLOOP;
           point.yCol _ m;
           colhead[m] _ [point:point, tag: y]};
        a[n][point.yCol] _ a[n][point.yCol] + yCoef

     END;
END;

EnterVer: PROC[i, j: PointPtr] =

BEGIN
   Constraint[i.x - j.x, eq];
   Partial[i, 1, 0]; Partial[j, -1, 0]
END;

EnterHor: PROC[i, j: PointPtr] =

BEGIN
   Constraint[i.y - j.y, eq];
   Partial[i, 0, 1]; Partial[j, 0, -1]
END;

EnterPara: PROC[i, j, k, l: PointPtr] =

BEGIN
   Constraint[(j.x - i.x)*(k.y - l.y) + (l.x - k.x)*(j.y - i.y), eq];
   Partial[i, l.y - k.y, k.x - l.x];
   Partial[j, k.y - l.y, l.x - k.x];
   Partial[k, i.y - j.y, j.x - i.x];
   Partial[l, j.y - i.y, i.x - j.x]
 END;

EnterCong: PROC[i, j, k, l: PointPtr, frame: Frame] =

BEGIN
    IF frame # [NIL, NIL, NIL]
       AND frame.xP#NIL AND frame.yP#NIL THEN -- relativized congruence

    {ux: REAL = frame.xP.x - frame.org.x;
     uy: REAL = frame.xP.y - frame.org.y;
     vx: REAL = frame.yP.x - frame.org.x;
     vy: REAL = frame.yP.y - frame.org.y;

     -- The following computations essentially multiply the coordinates of
     --   the constrained points by the inverse of the frame matrix,
     --    except that the division by the frame's determinant is left out.

     -- This still preserves the zeros of the constraint, and is MUCH easier to
     --    compute and differentiate. It also avoids division by zero if the 
     --    frame becomes singular. However, in that case
     --    it will consider the constraint as being satisfied.


     iju: REAL = vy*(j.x-i.x) - vx*(j.y - i.y); -- coords of (j-i) rel frame.
     ijv: REAL = -uy*(j.x-i.x) + ux*(j.y - i.y);
     klu: REAL = vy*(l.x-k.x) - vx*(l.y - k.y); -- coords of (l-k) rel frame
     klv: REAL = -uy*(l.x-k.x) + ux*(l.y - k.y);

     error: REAL = 0.5 * (iju*iju + ijv*ijv- klu*klu - klv*klv); -- error function

     dijx: REAL = vy*iju - uy*ijv; -- der. of error w.r.t. (j.x-i.x)
     dijy: REAL = -vx*iju + ux*ijv;
     dklx: REAL = -vy*klu + uy*klv;
     dkly: REAL = vx*klu - ux*klv;

     dvx: REAL = -iju*(j.y-i.y) + klu*(l.y - k.y); -- der. w.r.t. vx
     dvy: REAL = iju*(j.x-i.x) - klu*(l.x-k.x);
     dux: REAL = ijv*(j.y - i.y) - klv*(l.y - k.y); -- der. w.r.t. ux
     duy: REAL = -ijv*(j.x-i.x) + klv*(l.x-k.x);

     Constraint[error, eq];
     Partial[i, -dijx, -dijy];
     Partial[j, dijx, dijy];
     Partial[k, -dklx, -dkly];
     Partial[l, dklx, dkly];
     Partial[frame.xP, dux, duy];
     Partial[frame.yP, dvx, dvy];
     Partial[frame.org, -dux-dvx, -duy-dvy]} -- Phew!

    ELSE -- frame is orthogonal; absolute congruence

    {ijx: REAL = j.x - i.x;
     ijy: REAL = j.y - i.y;
     klx: REAL = l.x - k.x;
     kly: REAL = l.y - k.y;

     Constraint[0.5 * (ijx*ijx + ijy*ijy - klx*klx - kly*kly), eq];

     Partial[i, -ijx, -ijy];
     Partial[j, ijx, ijy];
     Partial[k, klx, kly];
     Partial[l, -klx, -kly]}

END;

EnterPerp: PROC[i, j, k, l: PointPtr, frame: Frame] =

BEGIN
    IF frame # [NIL, NIL,  NIL]
      AND frame.xP#NIL AND frame.yP#NIL THEN -- relativized perpendicularity

    {ux: REAL = frame.xP.x - frame.org.x; -- frame vectors.
     uy: REAL = frame.xP.y - frame.org.y;
     vx: REAL = frame.yP.x - frame.org.x;
     vy: REAL = frame.yP.y - frame.org.y;

     -- The following computations essentially multiply the coordinates of
     --   the constrained points by the inverse of the frame matrix,
     --    except that the division by the frame's determinant is left out.

     -- This still preserves the zeros of the constraint, and is MUCH easier to
     --    compute and differentiate. It also avoids division by zero if the 
     --    frame becomes singular. However, in that case
     --    it will consider the constraint as being satisfied.

     iju: REAL = vy*(j.x-i.x) - vx*(j.y - i.y); -- coords of (j-i) rel frame. 
     ijv: REAL = -uy*(j.x-i.x) + ux*(j.y - i.y);
     klu: REAL = vy*(l.x-k.x) - vx*(l.y - k.y); -- coords of (l-k) rel frame
     klv: REAL = -uy*(l.x-k.x) + ux*(l.y - k.y);

     error: REAL = iju*klu+ijv*klv; -- error function for this constraint

     dijx: REAL = vy*klu - uy*klv; -- der. of error w.r.t. (j.x-i.x)
     dijy: REAL = -vx*klu + ux*klv;
     dklx: REAL = iju*vy - ijv*uy;
     dkly: REAL = -iju*vx + ijv*ux;

     dvx: REAL = -(j.y-i.y)*klu - iju*(l.y-k.y); -- der. w.r.t. vx
     dvy: REAL = (j.x-i.x)*klu + iju*(l.x-k.x);
     dux: REAL = -(j.y-i.y)*klv - ijv*(l.y-k.y); -- der. w.r.t. ux
     duy: REAL = (j.x-i.x)*klv + ijv*(l.x-k.x);

     Constraint[error, eq];
     Partial[i, -dijx, -dijy];
     Partial[j, dijx, dijy];
     Partial[k, -dklx, -dkly];
     Partial[l, dklx, dkly];
     Partial[frame.xP, dux, duy];
     Partial[frame.yP, dvx, dvy];
     Partial[frame.org, -dux-dvx, -duy-dvy]} -- Phew!

    ELSE -- frame is orthogonal; absolute perpendicularity

    {Constraint[(j.x - i.x)*(l.x - k.x) + (j.y - i.y)*(l.y - k.y), eq];
     Partial[i, k.x - l.x, k.y - l.y];
     Partial[j, l.x - k.x, l.y - k.y];
     Partial[k, i.x - j.x, i.y - j.y];
     Partial[l, j.x - i.x, j.y - i.y]}

END;

EnterAt: PROC[pp: PointPtr, x, y: REAL, frame: Frame] =

BEGIN
    xT, yT: REAL;
    IF frame # [NIL, NIL,  NIL] THEN
       {mFrame _ Alg.GetFrameMatrix [frame, mFrame];
        [xT, yT] _ Alg.TransformPoint[x, y, mFrame]}
    ELSE
       {xT _ x; yT _ y}; 

    Constraint[xT - pp.x, eq]; 

    Partial[pp, -1, 0];
    IF frame # [NIL, NIL,  NIL] THEN -- enter derivatives with respect to frame

     {Partial[frame.org, 1, 0];
      IF frame.xP # NIL THEN
         {Partial[frame.xP, x, 0]; Partial[frame.org, -x, 0];
          IF frame.yP # NIL THEN
              {Partial[frame.yP, y, 0]; Partial[frame.org, -y, 0]}
          ELSE
              {Partial[frame.xP, 0, -y]; Partial[frame.org, 0, y]}}};

    Constraint[yT - pp.y, eq];

    Partial[pp, 0, -1];
    IF frame # [NIL, NIL,  NIL] THEN -- enter derivatives with respect to frame

     {Partial[frame.org, 0, 1];
      IF frame.xP # NIL THEN
         {Partial[frame.xP, 0, x]; Partial[frame.org, 0, -x];
          IF frame.yP # NIL THEN
              {Partial[frame.yP, 0, y]; Partial[frame.org, 0, -y]}
          ELSE
              {Partial[frame.xP, y, 0]; Partial[frame.org, -y, 0]}}}

    END;

EnterCcw: PROC[i, j, k: PointPtr, frame: Frame] =

BEGIN
   fudge: REAL = 500.0; -- Greg's mysterious fudge factor
   IF frame # [NIL, NIL,  NIL] THEN ERROR; 
   Constraint[(j.y-i.y)*(k.x-i.x) - (j.x-i.x)*(k.y-i.y) + fudge, geq];
   Partial[i, k.y-j.y, j.x-k.x];
   Partial[j, i.y-k.y, k.x-i.x];
   Partial[k, j.y-i.y, i.x-j.x]
END;

BEGIN

 p: PointPtr _ points;
 UNTIL p = NIL DO
   p.oldx _ Real.RoundLI[p.x];
   p.oldy _ Real.RoundLI[p.y];
   p.xCol _ 0; p.yCol _ 0;
   p _ p.link
  ENDLOOP;

END;

success _ TRUE; impossible _ FALSE;

BEGIN

 c: Stor.ConstrPtr _ constrs;
 UNTIL c = NIL DO
   WITH c SELECT FROM
           cc: Stor.HorPtr => EnterHor[cc.i, cc.j];
           cc: Stor.VerPtr => EnterVer[cc.i, cc.j];
           cc: Stor.ParaPtr => EnterPara[cc.i, cc.j, cc.k, cc.l];
           cc: Stor.CongPtr => EnterCong[cc.i, cc.j, cc.k, cc.l, cc.frame];
           cc: Stor.PerpPtr => EnterPerp[cc.i, cc.j, cc.k, cc.l, cc.frame];
           cc: Stor.AtPtr => EnterAt[cc.p, cc.x, cc.y, cc.frame];
           cc: Stor.CcwPtr => EnterCcw[cc.i, cc.j, cc.k, cc.frame];
        ENDCASE => ERROR;
   c _ c.link
  ENDLOOP

END;

IF success THEN RETURN;

BEGIN

 i, j, r, c: INT;
 pp: REAL;

 -- r is the row in which we are pivoting.
 -- c is the column in which we are pivoting.
 -- i and j are miscellaneous row and column indices respectively
 -- p is the reciprocal of the pivot element; also used as temp for swapping.

 FOR r IN [1..n] DO
   IF rowhead[r].tag = eq
     OR rowhead[r].tag = geq AND b[r] < -.01 THEN
     BEGIN

  pp _ 0;
  c _ 0;
  FOR j IN [1 .. m] DO
    IF ABS[a[r][j]] >= pp AND colhead[j].tag IN [x .. y] THEN
      {c _ j; pp _ ABS[a[r][c]]}
   ENDLOOP;

  IF c = 0 OR ABS[a[r][c]] < .001 THEN
    IF ABS[b[r]] > .01
      THEN {impossible _ TRUE; LOOP} -- this row is an unsatisfiable constraint
      ELSE LOOP; -- row is identically zero; constraint redundant.

  BEGIN
     -- adjust r to account for inequalities:
     kk: INTEGER _ r; -- save r
     FOR i: INT IN [1..n] DO
       IF i # kk AND rowhead[i].tag = geq AND a[i][1] >= 0.0
        AND a[i][c] < 0.0 AND b[i] - a[i][c] * b[r] / a[r][c] < 0.0 THEN r _ i
      ENDLOOP
  END;

  BEGIN
     temp: TableauLabel = colhead[c];
     colhead[c] _ rowhead[r];
     rowhead[r] _ temp
  END;

  pp _ 1.0 / a[r][c]; a[r][c] _ 1.0;
  -- divide everything in pivot row by the pivot element:
  b[r] _ b[r]*pp;
  FOR j IN [1..m] DO
    a[r][j] _ a[r][j] * pp
   ENDLOOP;
  -- subtract from other rows:
  FOR i IN [1..n] DO
    IF i # r THEN
      {FOR j IN [1..m] DO
         IF j # c THEN -- for each a[i,j] outside the pivot row and column
           a[i][j] _ a[i][j] - a[i][c] * a[r][j]; -- note that a[r,j] was already * pp.
        ENDLOOP;
       b[i] _ b[i] - a[i][c] * b[r]}
   ENDLOOP;
  -- Finally process pivot column:
  FOR i IN [1..n] DO
    IF i # r THEN a[i][c] _ -a[i][c] * pp
   ENDLOOP;

     END
  ENDLOOP;

END;

moved _ FALSE;

BEGIN

 p: Stor.PointPtr;
 FOR i: INT IN [1..n] DO
   IF rowhead[i].tag IN [x .. y]  THEN-- row was neither redundant nor unsatisfiable
     {p _ rowhead[i].point;
      IF rowhead[i].tag = x THEN
        {p.x _ p.x + b[i];  IF p.oldx # Real.RoundLI[p.x] THEN moved _ TRUE}
      ELSE
        {p.y _ p.y + b[i];  IF p.oldy # Real.RoundLI[p.y] THEN moved _ TRUE}};
 ENDLOOP;

END;

END;

   nit: INT _ 0;
   someMoved: BOOL;
   success _ FALSE; impossible _ FALSE;
   MarkMobilePoints[TRUE];
   [success, impossible, someMoved] _ NewtonRaphson[];
   WHILE nit < 20 AND NOT (success OR impossible) DO
     IF displayingMotions THEN DisplayMotions[points];
     nit _ nit + 1;
     [success, impossible, someMoved] _ NewtonRaphson[];
     success _ success OR NOT someMoved
   ENDLOOP;
   MarkMobilePoints[FALSE];
   RETURN

END;   

END.
X
HerculesSolverImpl.mesa (ex OldjunoSolverImpl)

Last Edited by: Stolfi, February 29, 1984 4:25:00 am PST

Was OldjunoSolverImpl
Coded July 1981 by Greg Nelson
Last Edited by: Gnelson, October 11, 1983 7:00 pm

This is the solver for an interactive graphics program that allows a user to specify the positions of points by declaring geometrical contraints that the positions are to satisfy. The program solves constraints (which can be viewed as simultaneous polynomial equations) by an n-dimensional version of Newton's method for finding the root of a differentiable function.

Each rows and column of the matrix represents a single scalar variable. If the row or column is labelled with the TableauLabel (p, x), then the scalar variable that the row or column represents is the x coordinate of the point p; similarly (p, y) means the y coordinate.

The label (*, eq) means the scalar variable is a "slack variable" associated with some constraint; the variable should have the value zero in any solution.

The label (*, geq) means the scalar variable should have a value greater than or equal to zero in any solution.

-- Prepare points for iteration: reset xCol and yCol  and save old coordinates:

-- Fill in tableau a[][] and constant column b[] by enumerating current constraints:

-- Now invert matrix:

     -- Set c so a[r,c] is largest of a[r,1], a[r, 2], a[r, 3], excluding
     -- columns in which we have already pivoted.

     -- We will pivot at a[r,c], if it exists and is not too small, and is legal.
     -- First check that it is not too small:

     -- Next we must check that pivoting at a[r, c] will not ruin any
     -- inequalities that are currently satisfied.  If there are any such inequalities,
     --    we change r to the index of one of them
     --   (namely, the first that would become unsatisfied),
     --    and pivot there instead.

     -- Exchange row and column labels.

     -- Do the pivoting.
     -- >>> Possible bug here: what if a[r][c]<0? Dividing by it would reverse the 
     --      sign of the inequality. Is that OK?

-- Now move the points, and check if any really moved:

�Ê˜š
œ/™/™9š
œ™Jš
œ™Jš
œ1™1——Jš
œð™ð—š
Ïk
˜
Jšœ1�œl�œ
˜¥
—šœ�˜Jšœ�œ9�œ˜a—š
�˜Jš�œ[˜`Jšœ�œ
˜Jšœ�œ˜JšÏbœ�œ
�œ�œ
˜(Jšœ�œ�œ
˜Jšœ	�œ˜Jšœ	�œ˜Jš	œ	�œ�œ�œ
�œ
˜)Jš
œ
�œ�œ�œ
�œ�œ
�œ
˜>Jšœ�œ�œ
˜Jšœ�œ
�œ	˜Jšœ�œ˜šœ�œ�œ˜9Jš
œ�™�Jš
œœ
™œ
Jš
œp™p—Jšœ
�œ�œ˜+Jšœ
�œ�œ˜+Jšœ�œ�œ
˜#Jšœ�œ�œ
˜)JšÏnœ�	œ�œ�œ�œ�œ
�œ�œ
˜DJšŸœ�	œ�œ�œ�œ�œ
�œ�œ
˜NJšœ	�œ�œ˜+Jš
œ�œ�œ
Ïc$œ 
œ %œ .œ˜ï
Jšœ�œ
 œ�œ Bœ
 #Ðcr œ Nœ˜�šœ
Ÿœ�œ�œ Á
˜ÿ
Jšœ
�œ�œ
˜—šœ
Ÿ
œ�œ�œ Ü
˜�Jšœ
�œ�œ
˜—š
Ÿœ�œ
�œ0�œ�œ˜hšœ
�˜šœ
Ÿœ�œ�œ ?˜fJšœ
�œ�œ�œ
�œ�œ
�œ�œ
�œ
˜v—šœ
Ÿœ�œ 3˜aJšœ
�œ�œ�œ
�œ�œ�œ$�œ;�œ
�œ
˜Æ
—šœ
Ÿ
œ�œ
�œ�œ !œ ¡ ‹
œ Ïr¡ Á
˜íšœ
�˜Jšœ�œ <˜LJšœ
Ÿœ�œ	�œ�œ %œx °
¡ �œ†�œ
˜Ä	š	œ
Ÿœ�œ
�œ œ’˜¶Jšœ
�œ�œ
˜—š	œ
Ÿœ�œ
�œ œÉ˜æJšœ
�œ�œ
˜—šœ
Ÿœ�œ
�œ œS œµ˜¹Jšœ
�œ/�œ
˜9—Jšœ
Ÿœ�œ
�œ œq œÏ�œ�œ
˜ƒJšœ
Ÿœ�œ�œ �
�œ"�œ�œ
�œ/�œ
�œ
˜ËJš
œ
Ÿœ�œ�œ œ'�œ�œ
˜_š
œ
Ÿ
œ�œ�œ Dœ 
˜uJš'œ
�œ œ�œ�œ
�œ/�œ�œ�œ
�œ�œ	�œ�œ�œ�œ�œ
�œ�œ�œ�œ
˜¥—šœ
Ÿœ�œ �œ Î
˜„šœ
�œ�œ
�œ�œ
�œ�œ
�œ�œ�˜UJšœ	�œ�œ
 œ!�œ�œ
�œ(�œ�œ
�œ�œ�œ€
˜Â—Jšœ�œ
˜
šœ�œ
�œ�œ�˜+Jšœ	�œ�œ
 œ!�œ�œ
�œ(�œ�œ
�œ�œ�œ€
˜Â—Jšœ�œ�œ
˜—šœ
Ÿœ�œ˜!Jšœ
�œF�œ
˜P—šœ
Ÿœ�œ˜!Jšœ
�œF�œ
˜P—šœ
Ÿ	œ�œ˜(Jšœ
�œÛ
�œ
˜å
—šœ
Ÿ	œ�œ&˜6šœ
�œ�œ
�œ�œ�œ	�œ
�œ
�œ
�œ
�œ
 ˜mJš>œ
�œ&�œ&�œ&�œ# Fœ @œ Gœ Kœ Iœ 4œ :œ�œ" œ�œ-�œ" œ�œ0�œ1 œ�œ "œ�œ �œ �œ�œ% œ
�œ,�œ& œ�œ© ˜œ
—šœ�œ
 +˜5Jš	œ�œ�œ�œ�œÃ
˜¦—Jšœ
�œ
˜—šœ
Ÿ	œ�œ&˜6šœ
�œ�œ
�œ�œ�œ�œ
�œ
�œ
�œ
�œ
 ˜sJš@œ
�œ œ
�œ&�œ&�œ# Fœ @œ Gœ Kœ Iœ 4œ :œ�œ" œ
�œ-�œ" œ�œ0�œ &œ�œ "œ�œ �œ�œ �œ# œ
�œ,�œ# œ�œ¨ ˜Ÿ
—šœ�œ
 1˜;Jš
ϊ
˜ä
—Jšœ
�œ
˜—šœ
Ÿœ�œ�œ˜8šœ
�œ
�œ�œ
�œ�œ�œ�œo�œY�œ
�œ�œ�œ�œ
 *˜ÔJšœ'�œ�œ
�œI�œ�œ
�œN�œF˜´—šœ=�œ
�œ�œ�œ�œ
 *˜ˆ
Jšœ'�œ�œ
�œI�œ�œ
�œN�œE˜³—Jšœ�œ
˜	—šœ
Ÿœ�œ#˜2Jšœ
�œ�œ
 "œ�œ
�œ�œ�œ�œ
�œ¬
�œ
˜š—Jšœ(¢œ¢œ™Pšœ
�˜Jš	œ�œ�œ
�œj�œ
˜›
—Jšœ
�œ
˜Jšœ¢œ¢œ$™UJšœ�œ�œ
˜$šœ
�˜Jšœ �œ�œ
�œ�œ�œ
�œÑ�œ�œ�˜À—Jšœ
�œ
˜Jšœ
�œ	�œ
�œ
˜Jš
œ™šœ
�˜šœ�œ�œ *œ
 -œ
 Aœ
 Mœ�œ�œ�œ�œ�œ�œ
�œ�˜óš
œ|™|Jšœ�œ�œ
�œ�œ
�œ�œ�œ
�œ�œ�œ
˜•
—š
œ€
™€
Jšœ�œ�œ
�œ�œ�œ
�œ�œ�œ�œ
 +œ�œ
�œ 1˜Ñ
—š
œ²™²Jš!œ�œ )œ	�œ 
œ�œ�œ
�œ�œ�œ�œ�œ�œ�œ'�œ
�œ�œ
˜”—š
œ(™(Jšœ�œ^�œ
˜j—š
œž
™ž
Jš9œ( 8œ�œ�œ�œ�œ œ�œ�œ�œ�œ�œ�œ�œ�œ
�œ�œ
 4œ	 œ' &œ�œ*�œ !œ�œ�œ�œ�œ�œ�œ
˜à——Jšœ�œ�œ
˜—Jšœ
�œ
˜Jš
œ7™7Jšœ	�œ
˜šœ
�˜Jš$œ�œ�œ
�œ�œ�œ�œ� /œ"�œ�œ�œ�œ	�œ�œ�œ�œ	�œ�œ
˜ï—Jšœ
�œ
˜—Jšœ
�œ
˜—Jš"œ	�œ�œ�œ�œ�œ=�œ
�œ
�œ
�œ
�œ�œ�œ}�œ
�œ�œ�œ�˜²—Jšœ
�œ˜——Jš�œ
˜J™—…—Hv]M