program IntervalMilneMethod;
{$APPTYPE CONSOLE}
uses System.SysUtils, Winapi.Windows, IntervalArithmetic32and64;

type intervalfunction = function (const T, Y : interval) : interval;

function Milne (const Tk1, Yk1, Tk2, Yk2,
                      Tk3, Yk3, Yk4,
                      DeltaT, DeltaY, H : interval;
                const F, PSI            : intervalfunction) : interval;
var H3H, H5, Tk1H3H, Yk1H3H, iv1, iv2 : interval;
begin
  iv1:=2*F(Tk1,Yk1)-F(Tk2,Yk2);
  iv1:=iv1+2*F(Tk3,Yk3);
  iv1:=4*H*iv1/3; // the term with Fk-1, Fk-2 and Fk-3
  H3H.a:=-3*H.b;
  H3H.b:=H.b;
  Tk1H3H:=Tk1+H3H; // the first argument for PSI
  Yk1H3H:=Yk1+H3H*F(DeltaT,DeltaY); // the second argument for PSI
  iv2:=PSI(Tk1H3H,Yk1H3H);
  iv2:=27*iv2-251*iv2;
  H5:=H*H;
  H5:=H5*H5;
  H5:=H*H5; // h^5
  iv2:=H5*iv2/720; // the term with PSI
  Result:=Yk4+iv1-iv2
end {Milne};

var i, k, m                       : Integer;
    hour, min, sec, msec          : Word;
    amh, real_h, w, width         : Extended;
    Tk, Yk, Tk1, Yk1, Tk2, Yk2,
    Tk3, Yk3, Tk4, Yk4, H, DeltaT,
    DeltaY                        : interval;
    hh                            : AnsiString;
    left, lib_name, right, s      : string;
    sf                            : Char;
    find                          : Boolean;
    FTY, PSITY                    : function (const T, Y : interval) : interval;
    ftxt                          : TextFile;
    time1, time2                  : TDateTime;
    ident                         : HMODULE;

begin
  Writeln;
  Writeln ('Solving the initial value problem y'' = f(t,y(t)), y(0) = y0');
  Writeln;
  Writeln ('for delta_t = {t in R : 0 <= t <= a}');
  Writeln ('and delta_y = {y in R : b_left <= y <= b_right}');
  Writeln;
  Writeln ('by an interval version of Milne''s method.');
  Writeln;
  Writeln ('This method requires a dynamic link library (DLL) with '
           +'definitions');
  Writeln ('of the following interval functions:');
  Writeln;
  Writeln ('F(T,Y)   - an interval extension of f(t,y(t)) with the header');
  Writeln ('             function FTY (const T, Y : interval) : interval;');
  Writeln ('PSI(T,Y) - an interval extension of psi(t,y(t)) = f^(3)(t,y(t)) = '
           +'y^(4)(t)');
  Writeln ('           with the header');
  Writeln ('             function PSITY (const T, Y : interval) : interval;');
  Writeln;
  repeat
    Write ('Do you have such a DLL (y/n)?: ');
    Readln (sf)
  until (sf='y') or (sf='n');
  if sf='n'
    then Exit;
  Writeln;
  repeat
    Writeln ('Enter the name of your DLL');
    Write ('(the extension .DLL will be added automatically): ');
    Readln (lib_name)
  until (lib_name<>'');
  lib_name:=lib_name+'.DLL';
  if not FileExists(lib_name)
    then begin
           Writeln;
           Writeln ('The file '+lib_name+' not found. The program will be '
                    +'closed.');
           Readln;
           Exit
         end;
  ident:=LoadLibrary(PWideChar(lib_name));
  if ident=0
    then begin
           Writeln;
           Writeln ('Your DLL has not been loaded. The program will be '
                    +'closed.');
           Readln;
           Exit
         end
    else begin
           Writeln;
           Writeln ('Your DLL has been loaded.')
         end;
  Writeln;
  find:=True;
  @FTY:=GetProcAddress(ident, 'FTY');
  if @FTY=nil
    then begin
           Writeln ('The function FTY has not been found in your DLL.');
           find:=false
         end
    else Writeln ('The function FTY has been found in your DLL.');
  if find
    then begin
           @PSITY:=GetProcAddress(ident,'PSITY');
           if @PSITY=nil
             then begin
                    Writeln ('The function PSITY has not been found in your '
                             +'DLL.');
                    find:=False
                  end
             else Writeln ('The function PSITY has been found in your DLL.')
         end;
  if not find
    then begin
           FreeLibrary (ident);
           Writeln;
           Writeln ('The program will be closed.');
           Readln;
           Exit
         end;
  Writeln;
  Writeln ('Enter data for your problem:');
  Writeln;
  Write ('delta_t = {t in R: 0 <= t <= a}                   a = ');
  Readln (w);
  amh:=w;
  Str (w:26, s);
  DeltaT.a:=0;
  DeltaT.b:=left_read(s);
  Write ('delta_y = (y in R: b_left <= y <= b_right)   b_left = ');
  Readln (w);
  Str (w:26, s);
  DeltaY.a:=left_read(s);
  Write ('                                            b_right = ');
  Readln (w);
  Str (w:26, s);
  DeltaY.b:=right_read(s);
  repeat
    Write ('Step size (>0 & <1):                              h = ');
    Readln (hh);
    Val (hh, real_h, m)
  until (m=0) and (real_h>0) and (real_h<1);
  Write ('Y0 = [left_y, right_y]                       left_y = ');
  Readln (w);
  Str (w:26, s);
  Yk4.a:=left_read(s);
  Write ('                                            right_y = ');
  Readln (w);
  Str (w:26, s);
  Yk4.b:=right_read(s);
  Write ('Y1 = [left_y, right_y]                       left_y = ');
  Readln (w);
  Str (w:26, s);
  Yk3.a:=left_read(s);
  Write ('                                            right_y = ');
  Readln (w);
  Str (w:26, s);
  Yk3.b:=right_read(s);
  Write ('Y2 = [left_y, right_y]                       left_y = ');
  Readln (w);
  Str (w:26, s);
  Yk2.a:=left_read(s);
  Write ('                                            right_y = ');
  Readln (w);
  Str (w:26, s);
  Yk2.b:=right_read(s);
  Write ('Y3 = [left_y, right_y]                       left_y = ');
  Readln (w);
  Str (w:26, s);
  Yk1.a:=left_read(s);
  Write ('                                            right_y = ');
  Readln (w);
  Str (w:26, s);
  Yk1.b:=right_read(s);
  find:=False;
  repeat
    Write ('maximum nunmber of steps (m*h <= a and m > 3)     m = ');
    Readln (m);
    if (m*real_h<=amh) and (m>3)
      then find:=True
  until find;
  repeat
    Write ('output every k steps (k <= m)           output_step = ');
    Readln (k)
  until (k>0) and (k<=m);
  Writeln;
  repeat
    Write ('Output (s - screen, f - file) : ');
    Readln (sf)
  until (sf='s') or (sf='f');
  if sf='f'
    then begin
           AssignFile (ftxt, 'IMM_RESULTS.txt');
           Rewrite (ftxt);
           Write (ftxt, 'SOLVING THE INITIAL VALUE PROBLEM BY AN INTERVAL ');
           Writeln (ftxt, 'VERSION OF MILNE''S METHOD');
           Writeln (ftxt, ' ');
           Writeln (ftxt, 'Step size:            h = ', hh);
           Writeln (ftxt, 'Number of steps:      m = ', m);
           Writeln (ftxt, ' ');
           Writeln ('Output to IMM_RESULTS.txt in current directory')
         end
    else Writeln ('Output to the screen');
  Writeln;
  H:=int_read(hh);
  Tk4:=0;     // t0 = tk-4
  Tk3:=Tk4+H; // t1 = tk-3
  Tk2:=Tk3+H; // t2 = tk-2
  Tk1:=Tk2+H; // t3 = tk-1
  if sf='s'
    then begin
           Writeln ('Initial data:');
           Writeln ('-------------');
           Writeln;
           Write ('     T = [')
         end
    else begin
           Writeln (ftxt, 'Initial data:');
           Writeln (ftxt, '-------------');
           Writeln (ftxt, ' ');
           Write (ftxt, '     T = [')
         end;
  iends_to_strings (Tk4, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Write ('     Y = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Write (ftxt, '     Y = [')
         end;
  iends_to_strings (Yk4, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Writeln;
           Write ('     T = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Writeln (ftxt, ' ');
           Write (ftxt, '     T = [')
         end;
  iends_to_strings (Tk3, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Write ('     Y = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Write (ftxt, '     Y = [')
         end;
  iends_to_strings (Yk3, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Writeln;
           Write ('     T = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Writeln (ftxt, ' ');
           Write (ftxt, '     T = [')
         end;
  iends_to_strings (Tk2, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Write ('     Y = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Write (ftxt, '     Y = [')
         end;
  iends_to_strings (Yk2, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Writeln;
           Write ('     T = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Writeln (ftxt, ' ');
           Write (ftxt, '     T = [')
         end;
  iends_to_strings (Tk1, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Write ('     Y = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Write (ftxt, '     Y = [')
         end;
  iends_to_strings (Yk1, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Writeln;
           Write ('DeltaT = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Writeln (ftxt, ' ');
           Write (ftxt, 'DeltaT = [')
         end;
  iends_to_strings (DeltaT, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Write ('DeltaY = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Write (ftxt, 'DeltaY = [')
         end;
  iends_to_strings (DeltaY, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Writeln;
           Writeln ('Solution:');
           Writeln ('---------');
           Writeln
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Writeln (ftxt, ' ');
           Writeln (ftxt, 'Solution:');
           Writeln (ftxt, '---------');
           Writeln (ftxt, ' ');
           Writeln ('Solving ...')
         end;
  time1:=Time;
  for i:=4 to m do
    begin
      Tk:=Tk1+H;
      Tk2:=Tk1-H;
      Tk3:=Tk2-H;
      Yk:=Milne(Tk1,Yk1,Tk2,Yk2,Tk3,Yk3,Yk4,DeltaT,DeltaY,H,FTY,PSITY);
      if i mod k = 0
        then begin
               iends_to_strings (Tk, left, right);
               if sf='s'
                 then begin
                        Write ('     T = [');
                        Writeln (left, ',', right, ']');
                        Write ('     Y = [')
                      end
                 else begin
                        Write (ftxt, '     T = [');
                        Writeln (ftxt, left, ',', right, ']');
                        Write (ftxt, '     Y = [')
                      end;
               iends_to_strings (Yk, left, right);
               if sf='s'
                 then Writeln (left, ',', right, ']')
                 else Writeln (ftxt, left, ',', right, ']');
               width:=int_width(Yk);
               if sf='s'
                 then Writeln ('         width =', width:11)
                 else Writeln (ftxt, '         width =', width:11);
               if sf='s'
                 then Writeln
                 else Writeln (ftxt, ' ')
             end;
      if i<m
        then begin
               Tk4:=Tk3;
               Yk4:=Yk3;
               Tk3:=Tk2;
               Yk3:=Yk2;
               Tk2:=Tk1;
               Yk2:=Yk1;
               Tk1:=Tk;
               Yk1:=Yk
             end
    end;
  FreeLibrary (ident);
  time2:=Time;
  DecodeTime (time2-time1, hour, min, sec, msec);
  if sf='s'
    then begin
           Writeln;
           Write ('Time of calculations: ', hour, ':');
           if min<10
             then Write ('0', min, ':')
             else Write (min, ':');
           if sec<10
             then Writeln ('0' , sec, '.', msec)
             else Writeln (sec, '.', msec)
         end
    else begin
           Writeln (ftxt, ' ');
           Write (ftxt, 'Time of calculations: ', hour, ':');
           if min<10
             then Write (ftxt, '0', min, ':')
             else Write (ftxt, min, ':');
           if sec<10
             then Writeln (ftxt, '0' , sec, '.', msec)
             else Writeln (ftxt, sec, '.', msec);
           CloseFile (ftxt);
         end;
  Writeln;
  Writeln ('End of calculations');
  Readln
end.