program Endversion;

{ f(x,s,t) = s + Arctan(xy') + pi }

use i_ari, mvi_ari, mv_ari;



function f_kappa(y0:interval):interval;
var pihalbe:interval;
begin
pihalbe:= 2*Arctan(1);
f_kappa:=pihalbe;
end;

function f_Konst(y0:interval):interval;
var pi:interval;
begin 
pi:=4*Arctan(1);
f_Konst:=pi + y0; 
end;

function f_Lipp(y0:interval;a:real):interval;
begin
f_Lipp:=a;
end;


function f(x,s,t:interval):interval;
var pi :interval;
begin
pi:=4*Arctan(1);
f:=s+Arctan(x*t)+pi;
end; 

function f_x(x,s,t:interval): interval;
begin
f_x:=t/(1 + sqr(x*t));
end;

function f_s(x,s,t:interval):interval;
begin
f_s:=1;
end;


function f_t(x,s,t:interval):interval;
begin
f_t:=x/(1 + sqr(x*t));
end;





procedure G(x:rvector;ix,q:ivector;var init:boolean;var slopes:ivector);
label 1;
var argminy,argmaxz        : rvector[lbound(q)..ubound(q)];
i,j                        : integer;
nenner,zaehler,links       : interval;
wesen, hilf                : interval;


begin{function}



i:=1;
while i<ubound(q)+1 do
begin{while}

{ neue Idee }
if ((diam(slopes[i])=0) and init) then {Zeile von slopes bleibt}
goto 1;

{ 1. Schritt: Bestimme argminy }

for j := lbound(q) to ubound(q) do argminy[j]:=ix[j].sup;


if i=lbound(q) then hilf:=q[lbound(q)] - argminy[lbound(q)+1] * 0.5  
else
if i=ubound(q) then hilf:=q[ubound(q)] - argminy[ubound(q)-1] * 0.5 
else
hilf:=q[i] - argminy[i-1] * 0.5 - argminy[i+1] *0.5 ;

if not (hilf.inf < 0) then {Einheitsvektor uebergeben}
                      slopes[i]:=0
else 
begin{hilf.inf < 0 => 2. Schritt}

{ 2. Schritt: Bestimme argmaxz }

for j := lbound(q) to ubound(q) do argmaxz[j]:=ix[j].inf;


if i=lbound(q) then hilf:=q[lbound(q)] - argmaxz[lbound(q)+1] * 0.5 
else
if i=ubound(q) then hilf:=q[ubound(q)] - argmaxz[ubound(q)-1] * 0.5 
else
hilf:=q[i] - argmaxz[i-1] * 0.5  - argmaxz[i+1] * 0.5 ;
if not (hilf.sup > 0) then slopes[i]:=-0.5 
else
begin{hilf.sup > 0 => 3. Schritt}

{ 3. Schritt }
if i=lbound(q) then hilf:=q[lbound(q)] - x[lbound(q)+1] * 0.5 
else
if i=ubound(q) then hilf:=q[ubound(q)] - x[ubound(q)-1] * 0.5
else
hilf:=q[i] - x[i-1] * 0.5 - x[i+1] * 0.5 ;

if hilf.inf>0 then begin { 3.1}
               zaehler:=hilf;

               if i=lbound(q) then 
               begin
               hilf:= ( argminy[lbound(q)+1] - x[lbound(q)+1] ) * 0.5;
               nenner:=hilf;
               links:=zaehler/nenner;
               wesen:=intval(links.inf,1);
               slopes[i]:=-0.5+wesen*0.5;
               end
               else
               if i=ubound(q) then 
               begin
               hilf:= ( argminy[ubound(q)-1] - x[ubound(q)-1] ) * 0.5;
               nenner:=hilf;
               links:=zaehler/nenner;
               wesen:=intval(links.inf,1);
               slopes[i]:=-0.5 + wesen*0.5;
               end
               else
               begin
               hilf:=(argminy[i+1] - x[i+1] + argminy[i-1] - x[i-1]) * 0.5;
               nenner:=hilf;
               links:=zaehler/nenner;
               wesen:=intval(links.inf,1);
               slopes[i]:=-0.5 + wesen*0.5;
               end;

               end {3.1}
else 
begin {3.2anfang}
if i=lbound(q) then hilf:=q[lbound(q)] - x[lbound(q)+1] * 0.5 
else
if i=ubound(q) then hilf:=q[ubound(q)] - x[ubound(q)-1] * 0.5  
else
hilf:=q[i] - x[i-1] * 0.5 - x[i+1] * 0.5 ;
if hilf.sup < 0 then begin {3.2eigentlich}
               zaehler:=hilf;

               if i=lbound(q) then 
               begin
               hilf:=(argmaxz[lbound(q)+1] - x[lbound(q)+1]) * 0.5;
               nenner:=hilf;
               links:=zaehler/nenner;
               wesen:=intval(links.inf,1);
               slopes[i]:=-wesen*0.5;
               end
               else
               if i=ubound(q) then 
               begin
               hilf:=(argmaxz[ubound(q)-1] - x[ubound(q)-1]) * 0.5;
               nenner:=hilf;
               links:=zaehler/nenner;
               wesen:=intval(links.inf,1);
               slopes[i]:= -wesen*0.5;
               end
               else
               begin
               hilf:=(argmaxz[i+1] - x[i+1] + argmaxz[i-1] - x[i-1]) * 0.5;
               nenner:=hilf;
               links:=zaehler/nenner;
               wesen:=intval(links.inf,1);
               slopes[i]:=-wesen*0.5;
               end;

     
               end {3.2eigentlich}
else
begin {3.3}
              wesen:=intval(0,1);
              slopes[i]:=-wesen*0.5;
              

end;{3.3}



end;{3.2anfang}

end;{3.Schritt}
end;{2.Schritt}
1: i:=i+1;
end;{while}
init:= true;

end;{function}

function igatri(g,b:ivector):ivector[lbound(b)..ubound(b)];

{ ------------------------------------------------------------------------- }

{ in diesem Programm kommen nur Tridiagonalmatrizen vor }

{ die LU Zerlegung reduziert sich auf }
{ for k=1 to n-1 do}
{ begin l[k+1,k]:=a[k+1,k]/a[k,k]; }
{       a[k+1,k+1]:= a[k+1,k+1] - l[k+1,k]*a[k,k+1]; }
{ end }

{ die Vorwaertszerlegung reduziert sich auf }
{ y[1]: b[1]; }
{ for i:= 2 to n do }
{ y[i]:=b[i] - l[i,i-1]*y[i-1] }

{ die Rueckwaertssubstitution reduziert sich auf }
{ x[n]:= y[n]/a[n,n]; }
{ for i:= n-1 downto 1 do }
{ x[i]:=(y[i]-a[i,i+1]*x[i+1])/a[i,i] }

{ man braucht also nur zwei Vektoren ltri und utri }

{ ------------------------------------------------------------------------- }

var ltri,utri,y,x : ivector[lbound(b)..ubound(b)];
    i,j,k         : integer;

begin
utri[lbound(b)]:= 1;
for k:= lbound(b) to ubound(b)-1 do 
         begin
              ltri[k+1]:= g[k+1]/utri[k];
              utri[k+1]:= 1 - ltri[k+1]*g[k];
         end;



    {Loese Ly=b}
     y[lbound(b)]:=b[lbound(b)];
     for i := 2 to ubound(b) do y[i]:=b[i]-ltri[i]*y[i-1];


     {Loese nun [R]x=y}

     x[ubound(b)]:= y[ubound(b)]/utri[ubound(b)];
     for i:= ubound(b)-1  downto lbound(b) do x[i]:=(y[i]-g[i]*x[i+1])/utri[i];

for i:= lbound(b) to ubound(b) do igatri[i]:=x[i];

end;
   
function Hmin(x:rvector;q:ivector):ivector[lbound(q)..ubound(q)];
var 
    i            : integer;
    hilf         : interval;

begin

for i:= lbound(q) to ubound(q) do
begin

if i=lbound(q) then hilf:= x[lbound(q)] -  x[lbound(q)+1]*0.5 + q[lbound(q)]
else
   if i=ubound(q) then hilf:= x[ubound(q)] -  x[ubound(q)-1]*0.5 + q[ubound(q)]
   else hilf:= x[i] - x[i-1]*0.5 - x[i+1]*0.5 + q[i];
               
            
if x[i] > hilf.sup then Hmin[i]:=hilf
else
   if x[i] < hilf.inf then Hmin[i]:=x[i]
   else Hmin[i]:=intval(hilf.inf,x[i]);


end;
end;

function durchmesser(v : ivector):real;
var i            : integer;
    bisher, hilf : real;
begin
bisher:= diam(v[lbound(v)]);
for i:= lbound(v) + 1 to ubound(v) do
begin
     hilf:=diam(v[i]);
     if hilf > bisher then bisher:=hilf;
end;
durchmesser:=bisher;
end;



function imax(a	: interval):interval;
begin
   if (a.sup < 0) or (a.sup=0) then imax:=0
   else if (a.inf > 0) or (a.inf=0) then imax:=a
   else imax:=intval(0,a.sup);
end;



function drin(a,b : interval):boolean;
begin
if ((a.inf < b.inf) or (a.sup > b.sup))  then drin:=false
else drin:=true;
end;



function max( a,b : real):real;
begin
if a<b then max:= b
else        max:=a
end;

function betrag(a:real):real;
begin 
if a> 0 then betrag:= a
else betrag :=-a
end;


function leer(a,b : interval):boolean;
begin
if (a.inf <= b.sup) and (b.inf <= a.sup) then leer := false
else leer:= true;
end;


procedure haupt(n:integer;y0:interval;var left:real;var a:real);
var

     
      q,zstart,zalt,zneu        : ivector[1..n];
      schritt                   : interval;
      i,j, bis                  : integer;
      hilf,h1,h2                : interval;
      c,unten,oben,Fstrich      : real;
      arg1,arg2,arg3,part       : interval;
      Fausw,fx,fs,ft            : interval;
      kappa,Konst,Lipp          : interval;
      slopes                    : ivector[1..n];
      mitte                     : rvector[1..n];
      newton, Hminvector        : ivector[1..n];
      rechteseite               : ivector[1..n];
      count, zaehler, zaehl     : integer;
      initialisiert             : boolean;
    
begin

kappa:=f_kappa(y0);
Konst:=f_Konst(y0);
Lipp:=f_Lipp(y0,a);





zneu:=intval(0,y0.sup);

     schritt:=a/intval(n+1);

     arg1:=intval(0,a);
     arg2:=intval(0,y0.sup);
     hilf:=Lipp*y0 + Konst*a;
     arg3:=intval(-hilf.sup,0);

     Fausw:=f(arg1,arg2,arg3);

     fx:=f_x(arg1,arg2,arg3);
     fs:=f_s(arg1,arg2,arg3);
     ft:=f_t(arg1,arg2,arg3);


     hilf:= abs(fx + fs * arg3 + ft * Fausw );
     Fstrich:=hilf.sup;

     h1:=sqr(schritt)*intval(0.5,1)*Fausw;
     h2:=schritt*sqr(schritt)*intval(-Fstrich,Fstrich)/2;
     for i:= 1 to n do q[i]:=(h1+h2)/2;
     q[1]:=q[1]-y0/2;


zaehler:=0;
slopes:=0;
initialisiert:=false;
repeat
  zaehler:= zaehler+1;
  zalt:=zneu;
     
  for i:= 1 to n do mitte[i]:=zalt[i].inf;
       
  G(mitte,zalt,q,initialisiert,slopes);

     count:=0;
     for i:= 1 to n do
          if not (diam(slopes[i])=0) then count:=count+1;
     writeln('slopes has ',count,' nondegenerate intervals');
 
  Hminvector:=Hmin(mitte,q);

  rechteseite:=igatri(slopes,Hminvector);
     
  newton:= mitte - rechteseite;

  for i:= 1 to n do zneu[i]:= zalt[i]**newton[i];


     if zaehler mod 10 = 0 then {neues [q]}
     begin
   
         bis:= 0;
         while (bis < n) and (zneu[bis+1].inf > 0) do bis:= bis + 1;


     if bis > 1 then 
     begin
      
     for i:= 1 to bis-1 do
     begin
     arg1:=i*schritt;
     arg2:=zneu[i];
     hilf:=Lipp*arg2 + Konst*(a-arg1);
     arg3:=intval(-hilf.sup,0);

     Fausw:=f(arg1,arg2,arg3); 

     if i> 1 then
     begin
     arg1:=intval((i-1)*<schritt.inf,(i+1)*>schritt.sup);
     arg2:=intval(zneu[i+1].inf,zneu[i-1].sup);
     hilf:=Lipp*zneu[i-1] + Konst*(a-(i-1)*schritt);
     arg3:=intval(-hilf.sup,0);
     end
     else
     begin
     arg1:=intval(0,2*>schritt.sup);
     arg2:=intval(zneu[2].inf,y0.sup);
     hilf:=Lipp*y0 + Konst*a;
     arg3:=intval(-hilf.sup,0);
     end;

     fx:=f_x(arg1,arg2,arg3);
     fs:=f_s(arg1,arg2,arg3);
     ft:=f_t(arg1,arg2,arg3);
     part:=f(arg1,arg2,arg3);

     hilf:= abs(fx + fs * arg3 + ft * part );
     Fstrich:=hilf.sup;


     h1:=sqr(schritt)*Fausw;
     h2:=schritt*sqr(schritt)*intval(-Fstrich,Fstrich)/3;
     q[i]:=(h1+h2)/2;
     end;
     q[1]:=q[1]-y0/2;

     arg1:=intval((bis-1)*<schritt.inf,a);
     arg2:=intval(0,zneu[bis-1].sup);
     hilf:=Lipp*zneu[bis-1] + Konst*(a-(bis-1)*schritt);
     arg3:=intval(-hilf.sup,0);

     Fausw:=f(arg1,arg2,arg3); 

     fx:=f_x(arg1,arg2,arg3);
     fs:=f_s(arg1,arg2,arg3);
     ft:=f_t(arg1,arg2,arg3);


     hilf:= abs(fx + fs * arg3 + ft * Fausw );
     Fstrich:=hilf.sup;


     h1:=sqr(schritt)*intval(0.5,1)*Fausw;
     h2:=schritt*sqr(schritt)*intval(-Fstrich,Fstrich)/2;
     for i:= bis to n do q[i]:=(h1+h2)/2;
     
     end
     else writeln('no improvement possible yet');


     end;

         
    
until zalt=zneu;


for i:= 1 to n do
             begin
               unten:=zneu[i].inf;
               oben:=zneu[i].sup;
               writeln('        [',unten:23:0:-1,',',oben:23:0:1,']');
             end;  
bis:= 0;
while (bis < n) and (zneu[bis+1].inf > 0) do bis:= bis + 1;
c:=bis*<schritt.inf;
if c > left then left:=c;
writeln('The free boundary c is greater than ',left:23:0:-1);
writeln('The free boundary c is less than ',a:23:0:1);
writeln(zaehler,' iteration steps have been done');
writeln('There were still ',count,' nondegenerate intervals');
                              
{ jetzt wollen wir mal sehen, ob wir a verbessern koennen }

bis:=n+1;
while (bis > 1) and (zneu[bis-1].sup =0) do bis:= bis - 1;
if bis < n+1 then a:=bis*>schritt.sup;
writeln('The free boundary c is less than ',a:23:0:1);           

end;


var n             : integer;
    y0,hilf,kappa : interval;
    a,sicherheit  : real;
    left          : real;

begin



hilf:=10;
y0:=1/hilf;
writeln('How many components do you want to consider?');
read(n);


kappa:=f_kappa(y0);
hilf:=sqrt(2*y0/kappa);
a:=hilf.sup;
left:=0;
repeat
sicherheit:=a;
haupt(n,y0,left,a);
until sicherheit=a;
         
end.