program vers;

use i_ari, mvi_ari, mv_ari;

type variablen = record
                 art: (w,z);
                 stelle: integer;
                 end;

     basis = dynamic array[*] of variablen;


procedure GJ(var A:imatrix; var q: ivector; pz,ps:integer);
var i,j : integer;
begin

q[pz]:=q[pz]/A[pz,ps];
for j := lb(A,2) to ub(A,2) do if j<> ps then A[pz,j]:=A[pz,j]/A[pz,ps];
A[pz,ps]:=1;


for i:= lb(A,1) to ub(A,1) do
if not (i = pz) then
              begin
                  for j:= lb(A,2) to ub(A,2) do
                  if j<>ps then A[i,j]:=A[i,j]-A[pz,j]*A[i,ps];
                  q[i]:=q[i] - q[pz]*A[i,ps];
              end;
for i:= lb(A,1) to ub(A,1) do if i<>pz then A[i,ps]:=0;

end;

function anf_pivotzeile(q:ivector):integer;

{ es wird angenommen, dass  nicht q>=0 gilt }


var
    stelle, i : integer;
    vorerst   : interval;

begin

stelle:=lb(q);
vorerst:=q[lb(q)];
for i:= lb(q)+1 to ub(q) do
if q[i].sup < vorerst.inf then begin
                               vorerst:= q[i];
                               stelle:=i;
                               end;
anf_pivotzeile:=stelle;
end;


function pivotzeile(A:imatrix;q:ivector;ps:integer):integer;
var
  loesbar       : boolean;
  i,stelle,t    : integer;
  test, vorerst : interval;

begin


loesbar:=false;
for i:= ub(q) downto lb(q) do
if A[i,ps].inf > 0 then
                   begin
                   loesbar:=true;
                   stelle:=i;
                   end;
if not loesbar then begin
                    writeln('Der Algorithmus kann keine Loesung finden');
                    pivotzeile:=0;
                    end
else
begin
vorerst:=q[stelle]/A[stelle,ps];
t:=stelle;
for i:= lb(q) to ub(q) do
if A[i,ps].sup > 0 then
                   begin
                   test:= q[i].inf/A[i,ps].sup;
                   if test.inf < vorerst.inf then
                           begin
                           vorerst:=test;
                           t:=i;
                           if A[t,ps].inf > 0 then stelle:=t;
                           end;
                   end;
if A[t,ps].inf <= 0 then pivotzeile:=stelle
else pivotzeile:=t;
end;
end;


procedure einlesen(var A:imatrix; var q : ivector);
var i,j : integer;
    M   : imatrix[lb(A,1)..ub(A,1),lb(A,1)..ub(A,1)];
begin
A:=0;
writeln('Geben Sie Ihre Matrix M zeilenweise ein!');
for i:= lb(M) to ub(M) do
for j:= lb(M) to ub(M) do
begin
read(M[i,j]);
A[i,ub(A,1)+j]:=-M[i,j];
end;
writeln('Rechte Seite q eingeben:');
for i:= lb(q) to ub(q) do
begin
     read(q[i]);
     A[i,i]:=1;
     A[i,ub(A,2)]:=-1
end;
writeln;


end;

procedure schreiben( A:imatrix; q : ivector);
var unten, oben : real;
    i,j         : integer;
begin

writeln('Als Tableau hat man derzeit:');

   for i:= lb(A,1) to ub(A,1) do
   begin
      for j:= lb(A,2) to ub(A,2) do
      begin
       {unten:=  A[i,j].inf;}
       {oben:=  A[i,j].sup;}
       {write('[',unten:2:1,',',oben:2:1,']');}
        write(A[i,j]);
      end;
      {unten:=  q[i].inf;}
      {oben:=   q[i].sup;}
      {writeln(' | [',unten:2:1,',',oben:2:1,']');}
       write(q[i]);
   end;
end;




procedure inverse(var a:imatrix);
var
    inv            :imatrix[lb(a,1)..ub(a,1),lb(a,1)..ub(a,1)];
    b,x,y          :ivector[lb(a,1)..ub(a,1)];
    i,j,k          :integer;
    hilf           :interval;



begin

for k:=lb(a,1) to ub(a,1)-1 do
begin

  for i:= k+1 to ub(a,1) do
  begin
      for j:= k+1  to ub(a,1) do
      begin
           a[i,j]:=a[i,j] - a[k,j] * a[i,k] /a[k,k];
      end;

      {jetzt wird die untere Dreiecksmatrix besetzt}
      a[i,k]:=a[i,k]/a[k,k];

  end;
end;




for k:=lb(a,1) to ub(a,1) do
begin{inverse}
for i:= lb(a,1) to ub(a,1) do if i=k then b[i]:=1 else b[i]:=0;

{Loese nun Ly=b}

for i:= lb(a,1) to ub(a,1) do
begin
    hilf:=0;
    for j:= lb(a,1) to i-1 do
    begin
        hilf:= hilf + a[i,j]*y[j];
    end;
    y[i]:=b[i]- hilf;
end;

{Loese nun Rx=y}

for i:= ub(a,1) downto lb(a,1) do
begin
     hilf:=0;
     for j:= i+1 to ub(a,1) do
     begin
         hilf:= hilf + a[i,j]*inv[j,k];
     end;
     inv[i,k]:= (y[i] - hilf)/a[i,i];
end;
end;{inverse}

for i:= lb(a,1) to ub(a,1) do
for j:= lb(a,1) to ub(a,1) do a[i,j]:=inv[i,j];

end{procedure};




function lemke(M:imatrix;q:ivector):ivector[lb(q)..ub(q)];

var
      A                  : imatrix[lb(q)..ub(q),lb(q)..2*ub(q)+1];
      B                  : basis[lb(q)..ub(q)];
      i,j,ps,pz          : integer;
      glob               : integer;
      loesung            : ivector[lb(q)..2*ub(q)];
      gegangen, jetzt    : variablen;
      vorsicht           : boolean;
      fertig             : boolean;


label 1;

{ Zentrale Idee des Algorithmus:                               }
{ Hat die Variable z_j die Basis verlassen,                    }
{ so ist im naechsten Schritt die Variable w_j                 }
{ die Pivotspalte und w_j kommt zur Basis                      }
{ Und zwar an die Stelle, die sich als Pivotzeile herausstellt }

{ In gegangen merken wir uns, die Variable, die gegangen ist   }
{ In jetzt, die die Pivotspalte wird und in die Basis kommt    }

begin

{ Initialisierung der Basis }

for i:= lb(q) to ub(q) do begin
                  B[i].art:=w;
                  B[i].stelle:=i;
                  end;

A:=0;
for i:=lb(q) to ub(q) do
for j:=lb(q) to ub(q) do A[i,ub(A,1)+j]:=-M[i,j];

for i:= lb(q) to ub(q) do
begin
     A[i,i]:=1;
     A[i,ub(A,2)]:=-1
end;


fertig := true;
for i:= lb(q) to ub(q) do if q[i].inf <0 then fertig:=false;
if fertig then goto 1;



{ erster pivotschritt }

ps:=2*ub(q)+1;
pz:= anf_pivotzeile(q);

glob:=pz;

GJ(A,q,pz,ps);


B[pz].art:= z;
B[pz].stelle:=0;

gegangen.art:=w;
gegangen.stelle:=pz;

{ z[pz] ist neue Pivotspalte }
{ w[pz] hat Basis verlassen }

while B[glob].stelle=0 do
begin

if gegangen.art=w then { Pivotspalte }
                  begin
                  ps:=ub(q)+gegangen.stelle;
                  jetzt.art:=z
                  end;
if gegangen.art=z then { Pivotspalte }
                  begin
                  ps:=gegangen.stelle;
                  jetzt.art:=w;
                  end;
jetzt.stelle:=gegangen.stelle;
{ Pivotspalte merken }

pz:=pivotzeile(A,q,ps);

GJ(A,q,pz,ps); { Basiswechsel }


if B[pz].art=z then gegangen.art:=z;
{ z[pz] hat Basis verlassen }

if B[pz].art=w then gegangen.art:=w;
{ w[pz] hat Basis verlassen }

gegangen.stelle:=B[pz].stelle;

{ Zur Basis kommt die Pivotspalte }

B[pz].art:=jetzt.art;
B[pz].stelle:=jetzt.stelle;



end;{while}


{ jetzt aufsammeln }

1: loesung:=0;
for i := lb(q) to ub(q) do
begin
if B[i].art=z then loesung[ub(q)+B[i].stelle]:=q[i];
if B[i].art=w then loesung[B[i].stelle]:=q[i];
end;


{ w - Mz =q koennte zwei Loesungen haben }
{ wovon eine nichtnegative Eintraege haben }
{ koennte. Die koennte in dem Aufgeblaehten q }
{ dann enthalten sein }

vorsicht:=false;
for i:= lb(q) to 2*ub(q) do if loesung[i].inf<0 then vorsicht:=true;
if vorsicht then writeln('the ILA fails !');

for i:= lb(q) to ub(q) do lemke[i]:=loesung[ub(b)+i];


end;




function vertex(A1,A2:imatrix;b:ivector):ivector[lb(b)..ub(b)];
var
    M,inv       : imatrix[lb(b)..ub(b),lb(b)..ub(b)];
    q,w,z       : ivector[lb(b)..ub(b)];
    i           : integer;
begin
inv:=A1;
inverse(inv);
M:= inv*A2;
q:=inv*b;
z:=lemke(M,q);
w:=q+M*z;
for i:=lb(b) to ub(b) do if z[i].sup =0 then vertex[i]:=w[i]
                         else vertex[i]:=-z[i];
end;


function huelle(A:imatrix;b:ivector):ivector[lb(b)..ub(b)];
var A1,A2                : imatrix[lb(b)..ub(b),lb(b)..ub(b)];
    b0,hull              : ivector[lb(b)..ub(b)];
    X                    : imatrix[1..2,1..4];
    i,j                  : integer;
    unten,oben           : real;

begin

{wir beginnen mit y1}

for j:= lb(b) to ub(b) do
begin
A1[1,j]:=A[1,j].inf;
A2[1,j]:=A[1,j].sup;
end;

for j:= lb(b) to ub(b) do
begin
A1[2,j]:=A[2,j].inf;
A2[2,j]:=A[2,j].sup;
end;


b0[1]:=b[1].sup;
b0[2]:=b[2].sup;


X[*,1]:=vertex(A1,A2,b0);




{jetzt y2}

for j:= lb(b) to ub(b) do
begin
A1[1,j]:=A[1,j].inf;
A2[1,j]:=A[1,j].sup;
end;

for j:= lb(b) to ub(b) do
begin
A2[2,j]:=A[2,j].inf;
A1[2,j]:=A[2,j].sup;
end;


b0[1]:=b[1].sup;
b0[2]:=b[2].inf;


X[*,2]:=vertex(A1,A2,b0);



{jetzt y3}

for j:= lb(b) to ub(b) do
begin
A2[1,j]:=A[1,j].inf;
A1[1,j]:=A[1,j].sup;
end;

for j:= lb(b) to ub(b) do
begin
A1[2,j]:=A[2,j].inf;
A2[2,j]:=A[2,j].sup;
end;


b0[1]:=b[1].inf;
b0[2]:=b[2].sup;


X[*,3]:=vertex(A1,A2,b0);


{jetzt y4}

for j:= lb(b) to ub(b) do
begin
A1[1,j]:=A[1,j].sup;
A2[1,j]:=A[1,j].inf;
end;

for j:= lb(b) to ub(b) do
begin
A1[2,j]:=A[2,j].sup;
A2[2,j]:=A[2,j].inf;
end;


b0[1]:=b[1].inf;
b0[2]:=b[2].inf;


X[*,4]:=vertex(A1,A2,b0);



hull[1].sup:=X[1,1].sup;
for j:= 2 to 4 do
if hull[1].sup < X[1,j].sup then hull[1].sup := X[1,j].sup;

hull[1].inf:=X[1,1].inf;
for j:= 2 to 4 do
if hull[1].inf > X[1,j].inf then hull[1].inf := X[1,j].inf;

hull[2].sup:=X[2,1].sup;
for j:= 2 to 4 do
if hull[2].sup < X[2,j].sup then hull[2].sup := X[2,j].sup;

hull[2].inf:=X[2,1].inf;
for j:= 2 to 4 do
if hull[2].inf > X[2,j].inf then hull[2].inf := X[2,j].inf;

huelle:=hull;



end;

function f(x:rvector):ivector[lb(x)..ub(x)];
begin
f[1]:= -x[1]*x[1] + x[2]*x[2] -1;
f[2]:= x[1]*x[1] - x[2];
end;

function Jac(x:ivector):imatrix[lb(x)..ub(x),lb(x)..ub(x)];
begin
Jac[1,1]:= -2*x[1];
Jac[1,2]:= 2*x[2];
Jac[2,1]:= 2*x[1];
Jac[2,2]:=-1;
end;

procedure haupt(n:integer);
var  mitte                  : rvector[1..n];
     start,xneu,xalt,fval   : ivector[1..n];
     jacobi                 : imatrix[1..n,1..n];
     oben,unten             : real;
     i,schritt              : integer;



begin

start:=intval(11,19)/10;
xneu:=start;
schritt:=0;
repeat
       schritt:=schritt+1;
       xalt:=xneu;
       mitte:=mid(xalt);
       jacobi:=Jac(xalt);
       fval:=f(mitte);
       xneu:= mitte - huelle(jacobi,fval);
       xneu:=xneu**xalt;

until xneu=xalt;

writeln('The solution is contained in:');
for i:= 1 to n do
begin
       unten:= xneu[i].inf;
       oben:=  xneu[i].sup;
       writeln('[',unten:23:0:-1,',',oben:23:0:1,']');
end;

writeln('We have used ',schritt,' Iteration steps');
end;


var n : integer;

begin



n:=2;
haupt(n);


end.