File : funs2.ads
with Ada.Text_IO, Ints, Reps, MiniFuns, MiniFuns.Ops, Taylors2, Zeros;
use Ada.Text_IO, Ints, Reps, MiniFuns;
pragma Elaborate (Ints,Reps,MiniFuns);
generic
PDeg: in Natural;
-------------------- begin_include scalars Scalar
type Scalar is private;
--- addition and multiplication
with function IsZero(S: Scalar) return Boolean is <>;
with procedure SetZero(S: out Scalar) is <>;
with procedure Neg(S: in out Scalar) is <>;
with procedure Neg(S1: in Scalar; S2: out Scalar) is <>;
with function "-"(S: Scalar) return Scalar is <>;
with procedure Add(S1: in Scalar; S2: in out Scalar) is <>;
with procedure Sum(S1,S2: in Scalar; S3: out Scalar) is <>;
with function "+"(S1,S2: Scalar) return Scalar is <>;
with procedure Sub(S1: in Scalar; S2: in out Scalar) is <>;
with procedure Diff(S1,S2: in Scalar; S3: out Scalar) is <>;
with function "-"(S1,S2: Scalar) return Scalar is <>;
with procedure Mult(I: in Integer; S: in out Scalar) is <>;
with procedure Prod(I: in Integer; S1: in Scalar; S2: out Scalar) is <>;
with function "*"(I: Integer; S: Scalar) return Scalar is <>;
with procedure MultAdd(I: Integer; S1: in Scalar; S2: in out Scalar) is <>;
with procedure Mult(S1: in Scalar; S2: in out Scalar) is <>;
with procedure Prod(S1,S2: in Scalar; S3: out Scalar) is <>;
with function "*"(S1,S2: Scalar) return Scalar is <>;
with procedure MultAdd(S1,S2: in Scalar; S3: in out Scalar) is <>;
with function "**"(S: Scalar; I: Integer) return Scalar is <>;
--- conversion and i/o
with function Scal(K: Integer) return Scalar is <>;
with function Scal(R: Rep) return Scalar is <>;
with procedure ShowType(S: in Scalar) is <>;
with procedure Show1(N: in String; S: in Scalar) is <>;
with procedure Show2(N: in String; S,S2: in Scalar) is <>;
with procedure Get(F: in File_Type; S: out Scalar; Decimal: in Boolean := False) is <>;
with procedure Put(F: in File_Type; S: in Scalar; Decimal: in Boolean := False) is <>;
-------------------- end_include scalars
-------------------- begin_include scalars-more Scalar
--- basic
with function IsNumeric(S: Scalar) return Boolean is <>;
with function IsSharp(S: Scalar) return Boolean is <>;
with function IsNumbers(S: Scalar) return Boolean is <>;
with procedure ModNumbers(S: in out Scalar) is <>;
with procedure ToBalls(S: in out Scalar) is <>;
with function Approx(S: Scalar) return Rep is <>;
--- sets
with function Center0(S: Scalar) return Boolean is <>;
with function Contains0(S: Scalar) return Boolean is <>;
with function Contains(S1,S2: Scalar) return Boolean is <>;
with function BallAt0(R: Rep) return Scalar is <>;
with function BallAt0(S: Scalar) return Scalar is <>;
with function DiskAt0(R: Rep) return Scalar is <>;
with function DiskAt0(S: Scalar) return Scalar is <>;
with function Disk(R1,R2: Rep) return Scalar is <>;
with function Disk(S: Scalar) return Scalar is <>;
with function Center(S: Scalar) return Scalar is <>;
with function CenterDisk(S: Scalar) return Scalar is <>;
with function ModCenter(S: Scalar) return Scalar is <>;
with function Union(S1,S2: Scalar) return Scalar is <>;
with function Intersection(S1,S2: Scalar) return Scalar is <>;
--- order
with function "<"(S1,S2: Scalar) return Boolean is <>;
with function "<="(S1,S2: Scalar) return Boolean is <>;
with function ">="(S1,S2: Scalar) return Boolean is <>;
with function ">"(S1,S2: Scalar) return Boolean is <>;
with function Min(S1,S2: Scalar) return Scalar is <>;
with function Max(S1,S2: Scalar) return Scalar is <>;
with function Inf(S: Scalar) return Rep is <>;
with function Sup(S: Scalar) return Rep is <>;
--- functions
with function Norm(S: Scalar) return Scalar is <>;
with function MaxNorm(S: Scalar) return Radius is <>;
with function Sqr(S: Scalar) return Scalar is <>;
with function "*"(R: Rep; S: Scalar) return Scalar is <>;
with function "/"(S: Scalar; R: Rep) return Scalar is <>;
with function Inv(S: Scalar) return Scalar is <>;
with function "/"(S1,S2: Scalar) return Scalar is <>;
with function Exp(S: Scalar) return Scalar is <>;
with function Log(S: Scalar) return Scalar is <>;
with function Sqrt(S: Scalar) return Scalar is <>;
-------------------- end_include scalars-more
with procedure UpMultAdd(R1,R2: in Rep; R3: in out Rep) is <>;
package Funs2 is -- functions of two variables
pragma Elaborate_Body;
package Taylors is new Taylors2 (PDeg => PDeg, Scalar => Scalar);
package MiniFuns_Ops is new MiniFuns.Ops (Scalar => Scalar);
subtype Taylor is Taylors.Taylor2;
subtype Affine is MiniFuns_Ops.Vec3; --- A*U+B*V+C (if U,V are the variables)
procedure MultAdd(S: in Scalar; P: in Affine; F: in out Taylor); -- add S*P
procedure Mult(P: in Affine; F: in out Taylor); -- P*F
procedure MultAdd(P: in Affine; F: in Taylor; G: in out Taylor); -- P*F
procedure Comp(U,V: in Affine; F: in Taylor; G: out Taylor); -- G := F(U,V)
type Fun is
record
A: Args := Zero_Args; -- defines U=(T+T0)*(T-T0)+B*(S-S0), V=A*(T+T0)*(T-T0)+(S-S0)
E: Boolean := True; -- even errors
P,Q: Taylor; -- represent P(U,V)+T*Q(U,V)
end record;
type Flop is array(1..2) of Fun; -- operators F -> L(1)*F+L(2)*Reflect(F)
function Rho1(F: Fun) return Radius; -- radius for 1st variable, U
function Rho2(F: Fun) return Radius; -- radius for 2nd variable, V
function Rho(F: Fun) return RadPair; -- domain radii
function Dom(F: Fun) return Domain; -- variables and radii
procedure SetRho(R: in RadPair; F: in out Fun); -- change radii
function IsZero(F: Fun) return Boolean; -- F=0
procedure SetZero(A: in Args; F: out Fun); -- F := 0, variables A
procedure SetZero(Template: in Fun; F: out Fun); -- F := 0, same variables as Template
procedure ToBalls(F: in out Fun); -- convert to balls
procedure Symm(F: in out Fun); -- extract symmetric part
function Symm(F: Fun) return Fun; -- extract symmetric part
procedure AddConst(C: in Scalar; F: in out Fun); -- F+C
procedure AddBall(C: in Scalar; F: in out Fun); -- F+BallAt0(C)
procedure AddT1(C: in Scalar; F: in out Fun); -- add C*T
procedure AddS1(C: in Scalar; F: in out Fun); -- add C*S
procedure AddX(C: in Scalar; F: in out Fun); -- add C*X
procedure AddY(C: in Scalar; F: in out Fun); -- add C*Y
function IsNumbers(F: Fun) return Boolean; -- known to be a polynomial
function ModNumbers(F: Fun) return Fun; -- part that may not be be polynomial
function Norm(R: RadPair; F: Fun) return Scalar; -- ||F||_R
function Norm(F: Fun) return Scalar; -- ||F||_Rho(F)
procedure Neg(F: in out Fun); -- -F
procedure Neg(F: in Fun; G: out Fun); -- -F
function "-"(F: Fun) return Fun; -- -F
procedure Add(F: in Fun; G: in out Fun); -- F+G
procedure Sum(F,G: in Fun; H: out Fun); -- F+G
function "+"(F,G: Fun) return Fun; -- F+G
procedure Sub(F: in Fun; G: in out Fun); -- F-G
procedure Diff(F,G: in Fun; H: out Fun); -- F-G
function "-"(F,G: Fun) return Fun; -- F-G
procedure Mult(S: in Scalar; F: in out Fun); -- S*F
procedure Prod(S: in Scalar; F: in Fun; G: out Fun); -- S*F
function "*"(S: Scalar; F: Fun) return Fun; -- S*F
procedure MultAdd(S: in Scalar; F: in Fun; G: in out Fun); -- add S*F
procedure Mult(F: in Fun; G: in out Fun); -- F*G
procedure Prod(F,G: in Fun; H: out Fun); -- F*G
function "*"(F,G: Fun) return Fun; -- F*G
function Sqr(F: Fun) return Fun; -- F*F
procedure MultAdd(F,G: in Fun; H: in out Fun); -- add F*G
procedure ChangeVar(A: in Args; F: in out Fun); -- change to variables A
procedure Scale(S: in Scalar; A: in Args; F: in out Fun); -- F(S.,S.)
procedure Scale(S: in Scalar; F: in out Fun); -- F(S.,S.)
--- derivatives ignore Q component
procedure MDer(R: in RadPair; F: in Fun; G: out Fun); -- generator of scalings
procedure DerTS(R: in RadPair; F: in Fun; F1,F2: out Fun); -- T,S derivatives
procedure Der(R: in RadPair; F: in Fun; F1,F2: out Fun); -- X,Y derivatives
procedure DerTS(R: in RadPair; F: in Fun; F1,F2,F11,F12,F22: out Fun); -- T,S derivatives
procedure Der(R: in RadPair; F: in Fun; F1,F2,F11,F12,F22: out Fun); -- X,Y derivatives
function ValTS(T,S: Scalar; F: Fun) return Scalar; -- T,S evaluation
function Val(X,Y: Scalar; F: Fun) return Scalar; -- X,Y evaluation
procedure ValDerTS(T,S: in Scalar; F: in Fun; D1,D2: out Scalar); -- T,S eval of derivatives
procedure ValDer(X,Y: in Scalar; F: in Fun; D1,D2: out Scalar); -- X,Y eval of derivatives
procedure ValDerTS(T,S: in Scalar; F: in Fun; D1,D2,D11,D12,D22: out Scalar); -- T,S eval of derivatives
procedure ValDer(X,Y: in Scalar; F: in Fun; D1,D2,D11,D12,D22: out Scalar); -- X,Y eval of derivatives
procedure SplitParity(F: in Fun; F0,F1: out Fun); -- split inot even and odd (in T)
procedure Reflect(F: in out Fun); -- change sign of F.Q
function Reflect(F: Fun) return Fun; -- change sign of F.Q
function "*"(L: Flop; F: Fun) return Fun; -- L*F
function "*"(L1,L2: Flop) return Flop; -- L*F
function "-"(L1,L2: Flop) return Flop; -- L1-L2
function Norm(R: RadPair; L: Flop) return Scalar; -- operator norm of L
function CompTSRadii(FT,FS: Fun; A: Args) return RadPair; -- radii needed for CompTS
function CompVWRadii(V: Fun; A: Args) return RadPair; -- radii needed for CompVW
function CompRadii(FX,FY: Fun; A: Args) return RadPair; -- radii needed for Comp
function CompXRadii(FX: Fun; A: Args) return RadPair; -- radii needed for CompX
function CompYRadii(FY: Fun; A: Args) return RadPair; -- radii needed for CompY
procedure CompTS(FT,FS: in Fun; A: in Args; FU,FV: out Fun); -- T,S composition for Args
procedure Comp(FX,FY: in Fun; A: in Args; FU,FV: out Fun); -- X,Y composition for Args
procedure CompTS(FT,FS,G: in Fun; H: out Fun); -- T,S composition
procedure Comp(FX,FY,G: in Fun; H: out Fun); -- X,Y composition: H := G(FX,FY)
procedure CompX(FX,G: in Fun; H: out Fun); -- H := G(FX,Y)
function CompX(FX,G: Fun) return Fun; -- G(FX,Y)
procedure CompY(FY,G: in Fun; H: out Fun); -- H := G(X,FY)
function CompY(FY,G: Fun) return Fun; -- G(X,FY)
function CompVW(V,F: Fun) return Fun; -- F(V,W), W(x,y)=-V(-y,-x)
function FF(Fx,Gy,V: Fun) return Fun; -- V -> Gy(x,V)+Fx(V,W), W=-SV
function DFF(Fxx,Fxy,Gyy,V: Fun) return Flop; -- derivative of FF
procedure Show1(N: in String; H: in Fun);
procedure Show2(N: in String; H1,H2: in Fun);
procedure Get(F: in File_Type; H: out Fun; Decimal: in Boolean := False);
procedure Put(F: in File_Type; H: in Fun; Decimal: in Boolean := False);
procedure Read(Name: in String; H: out Fun; Decimal: in Boolean := False);
procedure Write(Name: in String; H: in Fun; Decimal: in Boolean := False);
--- needs approximate solution in non-numeric mode
procedure RVal(G: in Fun; X,C: in Scalar; S: in out Scalar); -- solve G(x,s)=c
procedure LVal(F: in Fun; Y,C: in Scalar; S: in out Scalar); -- solve F(s,y)=c
procedure ZVal(F,G: in Fun; X,Y: Scalar; S: in out Scalar); -- solve F(x,z)+G(z,y)=0
procedure VVal(F,G: in Fun; X: Scalar; S: in out Scalar); -- solve G(x,s)+F(s,-s)=0
procedure YVal(G: in Fun; S: in out Scalar); -- solve Gx(s^2,s)=[Gy(s^2,s)]^2
ScalZero: Scalar renames Taylors.ScalZero;
ScalOne: Scalar renames Taylors.ScalOne;
ScalTwo: Scalar renames Taylors.ScalTwo;
ScalFour: Scalar renames Taylors.ScalFour;
Trunc: Boolean renames Taylors.Trunc;
TayZero: constant Taylor := Taylors.Zero_Taylor2;
private
subtype Components is Taylors.Components;
subtype TaylorMat is Taylors.Taylor2Mat;
function QuasiDeg2(D,I: Integer; C: Components) return Integer;
procedure Mult(P: in Affine; F: in out Taylor; Tmp: out Taylor);
procedure FillZero(D: in Natural; G: out TaylorMat);
procedure CheckArgs(F: in Fun; G: in out Fun; OK: out Boolean);
procedure FillZero(F: in out Fun);
procedure Mult(F: in Fun; G: in out Fun; Tmp: out Fun);
procedure DerTS(FA: in Args; DP: in TaylorMat; F1,F2: in out Fun);
procedure DerTS(FA: in Args; DP: in TaylorMat; F11,F12,F22: out Fun);
procedure Eval(V1,V2: in Fun; F: in Taylor; V3: out Fun);
pragma Inline (Rho1,Rho2,Rho);
package FunZeros is new Zeros (Number => Scalar);
use Taylors, MiniFuns_Ops, FunZeros;
T3: constant Integer := 3*Int(Trunc);
Eps1: constant Rep := 1.0E-05;
Eps2: constant Rep := 1.0E-16;
end Funs2;