Simplex Algorithm

[Albert Chan]

Another optimization, remove_eq_c(a,n) now also return pivots information.
simplex_le(a, ±n) use pivots to make less slack variables (smaller identity matrix)

Code:

pivot_index(a) :=
BEGIN
    LOCAL b := contains((a:=abs(a)), 1);
    IF b THEN RETURN b END
    IF (b:=max(a)) THEN b := contains(a,b) END
    RETURN b;
END;

remove_eq_c(a, n) :=
BEGIN
  LOCAL j,k, c := len(a[1]);
  LOCAL p := makelist(n);   /* max pivots */
  FOR j FROM n DOWNTO 1 DO
    p[j] := k := pivot_index(suppress(a[j],c));
    IF k==0 and a[j][c] THEN error(a) END
    a := k ? pivot((a[j]:=a[j]/a[j][k]),j,k) : delrows(a,j);
  END
  return a, remove(0,p);
END;

simplex_le(a, dir) :=
BEGIN
 LOCAL b, p := [];
 LOCAL n := abs(dir);
 IF n<len(a) THEN a,p := remove_eq_c(a,n) END
 b := dim(a);
 a := transpose(extend(a, identity(b(1))));
 a[0] := a[b[2]];
 a := transpose(delrows(a, b[2] .. b[2]+len(p)));
 n := sum(b)-2 - len(p);         /* variables */
 p := extend(p,range(b[2],n+1)); /* pivots */
 IF dir>0 THEN a[0] := -a[0] END /* maximize */
 a := simplex2(a,p,n,0,0);
 IF dir>0 THEN a[1] := -a[1] END /* undo flip */
 return a;
END;

Benchmark example, time reduced to 0.0050 s. For this example, no slack variables!

> m := [[3,2,1,10], [2,5,3,15], [-2,-3,-4,0]]
> simplex_le(m, -2)

\(\displaystyle [
\frac{-130}{7},[\frac{15}{7},0,\frac{25}{7}],\left(\begin{array}{cccc}
1 & \frac{1}{7} & 0 & \frac{15}{7} \\
0 & \frac{11}{7} & 1 & \frac{25}{7} \\
0 & \frac{25}{7} & 0 & \frac{130}{7}
\end{array}\right) ]\)

An Introduction to Linear Programming and Game Theory 3rd Edition, Paul R. Thie and Gerard E. Keough
Chapter 3.7 Redundant system, example 2. Again, no slack variables!

> m := [[1,2,0,1,20], [2,1,1,0,10], [-1,4,-2,3,40], [1,4,3,2,0]]
> simplex_le(m, -3)

\([35,[5,0,0,15],\left(\begin{array}{ccccc}
0 & \frac{3}{2} & \frac{-1}{2} & 1 & 15 \\
1 & \frac{1}{2} & \frac{1}{2} & 0 & 5 \\
0 & \frac{1}{2} & \frac{7}{2} & 0 & -35
\end{array}\right) ]\)

We have redundant "=" constraint, net of only 2 "=" constraints.

> remove_eq_c(m, 3)

\([\left(\begin{array}{ccccc}
0 & 1 & \frac{-1}{3} & \frac{2}{3} & 10 \\
1 & 0 & \frac{2}{3} & \frac{-1}{3} & 0 \\
0 & 0 & \frac{11}{3} & \frac{-1}{3} & -40
\end{array}\right) ,[2,1] ]\)

This simulate my previous version, not using pivot information.

> simplex_le(Ans[1], -inf)

\([35,[5,0,0,15,0,0],\left(\begin{array}{ccccccc}
0 & \frac{3}{2} & \frac{-1}{2} & 1 & \frac{3}{2} & 0 & 15 \\
1 & \frac{1}{2} & \frac{1}{2} & 0 & \frac{1}{2} & 1 & 5 \\
0 & \frac{1}{2} & \frac{7}{2} & 0 & \frac{1}{2} & 0 & -35
\end{array}\right) ]\)

Update: HP Prime had a weird behavior list /= k interpeted as list := list * (1/k)

> m := [3.0]
> m /= 3.0
> m[1] - 1      → −3.5527136788e−15

Thus, in remove_eq_c, to turn pivot as 1. (exactly), we do a[j] := a[j]/a[j][k]
I have made changes to ftneek code as well. Now pivot of 1. is exactly 1
sol(a,v) (and others) are reverted back to original, comparing pivot 1 instead of fuzzy 1.0