File : taylors2.ads

with Ada.Text_IO, Ints, Reps;
use Ada.Text_IO, Ints, Reps;

pragma Elaborate (Ints,Reps);

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 Scalar

-------------------- 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 Scalar

  with procedure UpMultAdd(R1,R2: in Rep; R3: in out Rep) is <>;

package Taylors2 is -- Taylor series in 2 variables, say U and V

  pragma Elaborate_Body;

  subtype Power is Integer range 0 .. PDeg;
  type Components is array(Power,Power) of Scalar;

  type Taylor2 is
    record
      R: RadPair;     --- domain radii
      C: Components;  --- coefficients and errors
    end record;

  type Taylor2Mat is array(Integer range <>,Integer range <>) of Taylor2;

  procedure SetRho(R: in RadPair; F: in out Taylor2);  -- change domain radius
  function IsZero(F: Taylor2) return Boolean;          -- F=0
  function Zero_Taylor2 return Taylor2;                -- 0
  procedure SetZero(R: in RadPair; F: out Taylor2);    -- F=0
  procedure ZapConst(F: in out Taylor2);               -- modulo constant
  procedure AddConst(C: in Scalar; F: in out Taylor2); -- F+C
  procedure UnitVector(I,J: in Natural; R: in RadPair; F: out Taylor2);     -- UnitVector(I,J)
  function IsNumbers(F: Taylor2) return Boolean;       -- known to be a polynomial
  function CenterDisk(F: Taylor2) return Taylor2;      -- part that is known to be polynomial
  procedure ModNumbers(F: in out Taylor2);             -- part that may not be be polynomial
  procedure ModNumbers(F: in Taylor2; G: out Taylor2); -- part that may not be be polynomial
  function ModNumbers(F: Taylor2) return Taylor2;      -- part that may not be be polynomial
  function Norm(R: RadPair; F: Taylor2) return Scalar; -- ||F||_R
  function Norm(F: Taylor2) return Scalar;             -- ||F||_Rho(F)
  procedure ToBalls(F: in out Taylor2);                -- erase polynomial info

  procedure Neg(F: in out Taylor2);                    -- -F
  procedure Neg(F: in Taylor2; G: out Taylor2);        -- -F
  function "-"(F: Taylor2) return Taylor2;             -- -F
  procedure Add(F: in Taylor2; G: in out Taylor2);     -- F+G
  procedure Sum(F,G: in Taylor2; H: out Taylor2);      -- F+G
  function "+"(F,G: Taylor2) return Taylor2;           -- F+G
  procedure Sub(F: in Taylor2; G: in out Taylor2);     -- F-G
  procedure Diff(F,G: in Taylor2; H: out Taylor2);     -- F-G
  function "-"(F,G: Taylor2) return Taylor2;           -- F-G
  procedure Mult(K: in Integer; F: in out Taylor2);                   -- K*F
  procedure Prod(K: in Integer; F: in Taylor2; G: out Taylor2);       -- K*F
  function "*"(K: Integer; F: Taylor2) return Taylor2;                -- K*F
  procedure MultAdd(K: in Integer; F: in Taylor2; G: in out Taylor2); -- add K*F to G
  procedure Mult(S: in Scalar; F: in out Taylor2);                    -- S*F
  procedure Prod(S: in Scalar; F: in Taylor2; G: out Taylor2);        -- S*F
  function "*"(S: Scalar; F: Taylor2) return Taylor2;                 -- S*F
  procedure MultAdd(S: in Scalar; F: in Taylor2; G: in out Taylor2);  -- add S*F

  procedure Mult(F: in Taylor2; G: in out Taylor2);                   -- F*G
  procedure Prod(F,G: in Taylor2; H: out Taylor2);                    -- F*G
  function "*"(F,G: Taylor2) return Taylor2;                          -- F*G
  procedure MultAdd(F,G: in Taylor2; H: in out Taylor2);              -- add F*G to H
  procedure Scale(S1,S2: in Scalar; F: in out Taylor2);               -- F(S1.S2.)
  function Val(S1,S2: Scalar; F: Taylor2) return Scalar;                           -- F(S1,S2)
  procedure Ders(R: in RadPair; F0: in Taylor2; F: out Taylor2);                   -- 0th derivative
  procedure Ders(R: in RadPair; F0: in Taylor2; F,Fu,Fv: out Taylor2);             -- derivatives
  procedure Ders(R: in RadPair; F0: in Taylor2; F,Fu,Fv,Fuu,Fuv,Fvv: out Taylor2); -- derivatives
  procedure Ders(D: in Natural; R: in RadPair; F: in Taylor2; G: out Taylor2Mat);  -- derivatives

  procedure Show1(N: in String; H: in Taylor2);
  procedure Show2(N: in String; H1,H2: in Taylor2);
  procedure Get(F: in File_Type; H: out Taylor2; Decimal: in Boolean := False);
  procedure Put(F: in File_Type; H: in Taylor2; Decimal: in Boolean := False);
  procedure Read(Name: in String; H: out Taylor2; Decimal: in Boolean := False);
  procedure Write(Name: in String; H: in Taylor2; Decimal: in Boolean := False);

  --- auxiliary
  function InDomain(R: RadPair; F: Taylor2) return Boolean;
  function QuasiDeg1(D: Integer; F: Taylor2) return Integer;
  function QuasiDeg(F: Taylor2) return Integer;

  ScalZero:    constant Scalar    := Scal(0);
  ScalOne:     constant Scalar    := Scal(1);
  ScalNegOne:  constant Scalar    := Scal(-1);
  ScalTwo:     constant Scalar    := Scal(2);
  ScalFour:    constant Scalar    := Scal(4);
  Trunc:       constant Boolean   := IsNumeric(ScalZero);

private

  procedure Collect(R1,R2: in Radius; DF: in Natural; F: in out RepMatrix; G: in out Components);
  procedure MultAdd(DF,DG: in Natural; F,G: in Taylor2; H: in out Taylor2);
  procedure PolyDerU(F: in Taylor2; G: out Taylor2);
  procedure PolyDerV(F: in Taylor2; G: out Taylor2);
  function DerFac1(R1,R2: Radius) return Scalar;
  function DerFac2(R1,R2: Radius) return Scalar;

  T3:    constant Integer   := 3*Int(Trunc);
  EZero: constant RepMatrix := ZeroRepMat(2*PDeg);
  Tay_Zero: Taylor2;

end Taylors2;