unit UInterv;
{ Interval Arithmetic Unit v0.5                          }
{ Delphi 2.0 Source                                      }
{ Made by Joao Paulo Schwarz Schuler                     }
{ http://www.schulers.com/jpss/pascal                    }
{ you should see the exmaples for help on the above site }
{
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Library General Public
License as published by the Free Software Foundation.

This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
Library General Public License for more details.
}

interface

type TIFloat = extended;  // Interval Float

type TInterval = array[1..2] of TIFloat;

const csIZero:TInterval = (0,0);
function ISoma(A,B:TInterval):TInterval;        // A+B
function INeg(A:TInterval):TInterval;           //  -A
function ISub(A,B:TInterval):TInterval;         // A-B
function IMul(A,B:TInterval):TInterval;         // A*B
function IInv(A:TInterval):TInterval;           // 1/A
function IDiv(A,B:TInterval):TInterval;         // A/B
function IUni(A,B:TInterval):TInterval;         // A U B
function IInter(A,B:TInterval):TInterval;       // A  B
function IDist(A,B:TInterval):TIFloat;
function IMed(A:TInterval):TIFloat;             // Medium
function FloatToInterval(F:TIFloat):TInterval;
function IntervalToStr(F:TInterval):String;
function IMod(A:TInterval):TIFloat;             // |A|
function IDiam(A:TInterval):TIFloat;            // Diameter
function ISet(X,Y:TIFloat):TInterval;           {define a interval}
function IIn(X:TIFloat;A:TInterval):boolean;   // In

type T1FFunction = function (X:TIFloat):TIFloat;
{ Unary Float Function }

function IApplyFn(F:T1FFunction;X:TInterval):TInterval;

type TIntervalNewtonSolution = object
                                 private
                                 F:T1FFunction; // function
                                 D:T1FFunction; // f'
                                 M,X:TInterval;

                                 public
                                 constructor Init(PF,PD:T1FFunction;PX:TInterval);
                                 { PF:f; PD: f' ; PN: Initial Interval}

                                 function RunStep:TInterval;
                                 function GetLast:TInterval;

                                 destructor Done;
                               end;

implementation
{$O-}

uses SysUtils; // it's a Delphi library

function IApplyFn(F:T1FFunction;X:TInterval):TInterval;
begin
IApplyFn:=ISet(F(X[1]),F(X[2]));
end;

constructor TIntervalNewtonSolution.Init(PF,PD:T1FFunction;PX:TInterval);
                                 { PF:f; PD: f' ; PX: Initial Interval}
begin
F:=PF;
D:=PD;
X:=PX;
end;

function TIntervalNewtonSolution.RunStep:TInterval;
var MX:TIFloat;
    IMX,NU,DE:TInterval;
begin
MX:=IMed(X);
IMX:=FloatToInterval(MX);
NU:=FloatToInterval(F(MX));
DE:= IApplyFn(D,X);

X:=IInter(X, ISub(IMX,IDiv( NU ,DE  ) ));
RunStep:=X;
end;

function TIntervalNewtonSolution.GetLast:TInterval;
begin
GetLast:=X;
end;

destructor TIntervalNewtonSolution.Done;
begin
end;

function MaxE(A,B:TIFloat):TIFloat;
begin
if B>A
   then MaxE:=B
   else MaxE:=A
end;

function MinE(A,B:TIFloat):TIFloat;
begin
if B<A
   then MinE:=B
   else MinE:=A
end;

procedure ExchangeE(var A,B:TIFloat);
var AUX:TIFloat;
begin
A:=AUX;
A:=B;
B:=AUX;
end;

procedure Ordena(var A,B:TIFloat);
begin
if B<A
   then ExchangeE(A,B);
end;

function IGetMinMax(X1,X2,X3,X4:TIFloat):TInterval;
begin
Ordena(X1,X2);Ordena(X3,X4);
Ordena(X1,X3);Ordena(X2,X4);
IGetMinMax[1]:=X1;
IGetMinMax[2]:=X4;
end;

function ISoma(A,B:TInterval):TInterval;   // A+B
begin
ISoma[1]:=A[1]+B[1];
ISoma[2]:=A[2]+B[2];
end;

function INeg(A:TInterval):TInterval;  //  -A
begin
INeg[1]:=-A[2];
INeg[2]:=-A[1];
end;

function ISub(A,B:TInterval):TInterval;  // A-B
begin
ISub:=ISoma(A,INeg(B));
end;

function IMul(A,B:TInterval):TInterval;  // A*B
begin
IMul:=IGetMinMax(A[1]*B[1],A[1]*B[2],A[2]*B[1],A[2]*B[2]);
end;

function IInv(A:TInterval):TInterval;  // 1/A
begin
IInv[1]:=1/A[2];
IInv[2]:=1/A[1];
end;

function IDiv(A,B:TInterval):TInterval;  // A/B
begin
IDiv:=IMul(A,IInv(B));
end;

function IUni(A,B:TInterval):TInterval;  // A U B
begin
IUni[1]:=MinE(A[1],B[1]);
IUni[2]:=MaxE(A[2],B[2]);
end;

function IInter(A,B:TInterval):TInterval;  // A  B
var X,Y:TIFloat;
begin
X:=MaxE(A[1],B[1]);
Y:=MinE(A[2],B[2]);
if Y<X
   then IInter[1]:=Y
   else IInter[1]:=X;
IInter[2]:=Y;
end; // of procedure

function IDist(A,B:TInterval):TIFloat;
begin
IDist:=MaxE(Abs(A[1]-B[1]),Abs(A[2]-B[2]));
end;

function IMed(A:TInterval):TIFloat;
begin
IMed:=(A[1]+A[2])/2;
end;

function FloatToInterval(F:TIFloat):TInterval;
begin
FloatToInterval[1]:=F;
FloatToInterval[2]:=F;
end;

function IMod(A:TInterval):TIFloat;   // |A|
begin
IMod:=IDist(A,FloatToInterval(0));
end;

function IDiam(A:TInterval):TIFloat;   // Diameter
begin
IDiam:=Abs(A[1]-A[2]);
end;

function ISet(X,Y:TIFloat):TInterval;
{define a interval}
begin
ISet[1]:=X;
ISet[2]:=Y;
end;

function IntervalToStr(F:TInterval):String;
begin
IntervalToStr:='['+FloatToStr(F[1])+';'+FloatToStr(F[2])+']';
end;

function IIn(X:TIFloat;A:TInterval):boolean;   // In
begin
IIn:=((A[1]<=X) and (X<=A[2]))
end;

end.