with Ada.Text_IO, Ints, Logicals, Scalars, Rationals, Flts;
use Ada.Text_IO, Ints, Logicals, Scalars, Rationals, Flts;
with Recycling, Matrices, Vectors, Polynoms, Polynoms.More;

pragma Elaborate_All (Ada.Text_IO,Logicals,Scalars,Rationals,Flts);
pragma Elaborate_All (Recycling,Matrices,Polynoms.More);

generic

  PDeg: in Natural;  --- max degree

  -------------------- begin_include scalars Scalar

  type Scalar is private;

  --- basic
  with function Info(Dummy: Scalar) return Scalar_Info is <>;
  with function IsSharp(S: Scalar) return Boolean is <>;
  with function IsZero(S: Scalar) return Boolean is <>;
  with function "="(S1,S2: Scalar) return Boolean is <>;
  with procedure SetZero(S: in out Scalar) is <>;
  with procedure Copy(S1: in Scalar; S2: in out Scalar) is <>;
  with procedure Swap(S1,S2: in out Scalar) is <>;

  --- sets
  with function Center0(S: Scalar) return Boolean is <>;
  with function Contains0(S: Scalar) return Logical is <>;
  with function Contains(S1,S2: Scalar) return Logical is <>;
  with procedure BallAt0(R: in Flt; S: in out Scalar) is <>;
  with function BallAt0(R: Flt) return Scalar is <>;
  with procedure ToErr(S: in out Scalar) is <>;
  with procedure ToErr(S1: in Scalar; S2: in out Scalar) is <>;
  with function ToErr(S: Scalar) return Scalar is <>;
  with procedure Center(S: in out Scalar) is <>;
  with procedure Center(S1: in Scalar; S2: in out Scalar) is <>;
  with function Center(S: Scalar) return Scalar is <>;
  with procedure ModCenter(S: in out Scalar) is <>;
  with procedure ModCenter(S1: in Scalar; S2: in out Scalar) is <>;
  with function ModCenter(S: Scalar) return Scalar is <>;
  with procedure ErrMult(R: in Radius; S: in out Scalar) is <>;
  with procedure Union(S1: in Scalar; S2: in out Scalar) is <>;
  with function Union(S1,S2: Scalar) return Scalar is <>;
  with procedure Intersection(S1: in Scalar; S2: in out Scalar; Empty: out Logical) is <>;

  --- order
  with function Sign(S: Scalar) return Integer is <>;
  with function Compare(S1,S2: Scalar) return Integer 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 ">="(S1,S2: Scalar) return Boolean is <>;
  with procedure Min(S1: in Scalar; S2: in out Scalar) is <>;
  with procedure Min(S1,S2: in Scalar; S3: in out Scalar) is <>;
  with function Min(S1,S2: Scalar) return Scalar is <>;
  with procedure Max(S1: in Scalar; S2: in out Scalar) is <>;
  with procedure Max(S1,S2: in Scalar; S3: in out Scalar) is <>;
  with function Max(S1,S2: Scalar) return Scalar is <>;
  with function Inf(S: Scalar) return Flt is <>; -- inf of real part
  with function Sup(S: Scalar) return Flt is <>; -- sup of real part

  --- arithmetic
  with procedure Neg(S: in out Scalar) is <>;
  with procedure Neg(S1: in Scalar; S2: in out Scalar) is <>;
  with function "-"(S: Scalar) return Scalar is <>;
  with procedure Add(I: in Integer; S: in out Scalar) is <>;
  with procedure Add(S1: in Scalar; S2: in out Scalar) is <>;
  with procedure Sum(S1,S2: in Scalar; S3: in 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: in out Scalar) is <>;
  with function "-"(S1,S2: Scalar) return Scalar is <>;
  with procedure Mult(R: in Flt; S,Tmp: in out Scalar) is <>;
  with procedure Mult(R: in Flt; S: in out Scalar) is <>;
  with procedure Prod(R: in Flt; S1: in Scalar; S2: in out Scalar) is <>;
  with function "*"(R: Flt; S: Scalar) return Scalar is <>;
  with procedure AddProd(R: in Flt; S1: in Scalar; S2,Tmp: in out Scalar) is <>;
  with procedure AddProd(R: in Flt; S1: in Scalar; S2: in out Scalar) is <>;
  with procedure Mult(Q: in Rational; S: in out Scalar) is <>;
  with procedure Mult(S1: in Scalar; S2: in out Scalar) is <>;
  with procedure Prod(S1,S2: in Scalar; S3: in out Scalar) is <>;
  with function "*"(S1,S2: Scalar) return Scalar is <>;
  with procedure AddProd(S1,S2: in Scalar; S3,Tmp: in out Scalar) is <>;
  with procedure AddProd(S1,S2: in Scalar; S3: in out Scalar) is <>;
  with procedure SubProd(S1,S2: in Scalar; S3,Tmp: in out Scalar) is <>;
  with procedure Div(R: in Flt; S,Tmp: in out Scalar) is <>;
  with procedure Div(R: in Flt; S: in out Scalar) is <>;
  with procedure Quot(S1: in Scalar; R: in Flt; S2: in out Scalar) is <>;
  with function "/"(S: Scalar; R: Flt) return Scalar is <>;
  with procedure Div(S1: in Scalar; S2: in out Scalar) is <>;
  with procedure Quot(S1,S2: in Scalar; S3: in out Scalar) is <>;
  with function "/"(S1,S2: Scalar) return Scalar is <>;
  with procedure Inv(S1: in Scalar; S2: in out Scalar) is <>;
  with function Inv(S: Scalar) return Scalar is <>;
  with function "**"(S: Scalar; I: Integer) return Scalar is <>;

  --- fun
  with function IsReal(S: Scalar) return Boolean is <>;
  with procedure Adjoint(S: in out Scalar) is <>;
  with procedure Adjoint(S1: in Scalar; S2: in out Scalar) is <>;
  with function Adjoint(S: Scalar) return Scalar is <>;
  with procedure Real_Part(S: in out Scalar) is <>;
  with procedure Imag_Part(S: in out Scalar) is <>;
  with procedure Eval0(S: in out Scalar) is <>;
  with function "abs"(S: Scalar) return Scalar is <>;
  with procedure Norm(S1: in Scalar; S2: in out Scalar) is <>;
  with function Norm(S: Scalar) return Scalar is <>;
  with function MaxNorm(S: Scalar) return Radius is <>;
  with procedure Sqr(S: in out Scalar) is <>;
  with function Sqr(S: Scalar) return Scalar is <>;
  with procedure Sqr(S1: in Scalar; S2: in out Scalar) is <>;
  with procedure Sqrt(S1: in Scalar; S2: in out Scalar) is <>;
  with function Sqrt(S: Scalar) return Scalar is <>;
  with procedure Root(K: in Positive; S1: in Scalar; S2: in out Scalar) is <>;
  with function Root(K: Positive; S: Scalar) return Scalar is <>;
  with procedure Exp(S1: in Scalar; S2: in out Scalar) is <>;
  with function Exp(S: Scalar) return Scalar is <>;
  with procedure ArcCos(S1: in Scalar; S2: in out Scalar) is <>;
  with function ArcCos(S: Scalar) return Scalar is <>;
  with procedure ArcSin(S1: in Scalar; S2: in out Scalar) is <>;
  with function ArcSin(S: Scalar) return Scalar is <>;
  with procedure Cos(S1: in Scalar; S2: in out Scalar) is <>;
  with function Cos(S: Scalar) return Scalar is <>;
  with procedure Sin(S1: in Scalar; S2: in out Scalar) is <>;
  with function Sin(S: Scalar) return Scalar is <>;
  with procedure Log(S1: in Scalar; S2: in out Scalar) is <>;
  with function Log(S: Scalar) return Scalar is <>;
  with function Pi return Scalar is <>;
  with procedure Simple_Random(S: in out Scalar) is <>;
  with procedure QPower(Q: in Rational; S1: in Scalar; S2: in out Scalar) is <>;
  with procedure Roots2(B,C: in Scalar; U1,U2,V: in out Scalar) is <>;
  with procedure Roots3(B,C,D: in Scalar; U0,U1,U2,V: in out Scalar) is <>;

  --- conversion and i/o
  with function Approx(S: Scalar) return Flt is <>;
  with procedure Assign(K: in Integer; S: in out Scalar) is <>;
  with procedure Assign(Q: in Rational; S: in out Scalar) is <>;
  with procedure Assign(R: in Flt; S: in out Scalar) is <>;
  with function Scal(I: Integer) return Scalar is <>;
  with function Scal(Q: Rational) return Scalar is <>;
  with function Scal(R: Flt) return Scalar is <>;
  with procedure Enclose(R1,R2: in Flt; S: in out Scalar) is <>;
  with function Enclose(R1,R2: Flt) return Scalar is <>;
  with procedure Show1(N: in String; S: in Scalar; Whatever: in Boolean := True) is <>;
  with procedure Show2(N: in String; S1,S2: in Scalar; Whatever: in Boolean := True) is <>;
  with procedure Put(F: in File_Type; S: in Scalar; Decimal: in Boolean := False) is <>;
  with procedure Get(F: in File_Type; S: in out Scalar; Decimal: in Boolean := False) is <>;

  --- misc
  with function Get_Precision(R: Scalar) return Positive is <>;
  with procedure Proper_Rounding is <>;
  with procedure Proper_Rounding(Dummy: in Scalar) is <>;

  -------------------- end_include scalars

  with package SM is new Matrices (Scalar => Scalar);
  with package SV is new Vectors (Component => Scalar);
  with package SP is new Polynoms (Scalar => Scalar, SV => SV);

package Taylors1 is  --- Taylor series with non-constant coefficients
                     --- see the WARNING below

  pragma Elaborate_Body;

  package SP_More is new SP.More;

  subtype Power is Integer range 0 .. PDeg;
  subtype Poly1 is SP.Polynom1(Power);

  Rad:               Radius  := One;    -- default domain
  Psi:               Radius  := 1.0E-3; -- default relative precision for Sup,Inf
  Use_FBalls:        Boolean := False;  -- use only when Scalar = MPFR.Floats.Balls.Ball
  Check_Consistency: Boolean := False;  -- optional checks in many procedures

  FConst: constant Positive := 16384;   -- needs to be larger than PDeg+1

  type Taylor1 is
    record
      C: Poly1;              -- coefficient C(I) with I<F have no general errors
      F: Natural := FConst;  -- this particular value means that P is constant
      R: Radius  := Rad;     -- domain of analyticity is |z|<Rad
    end record;

  --- WARNING: changing P.C directly can make P inconsistent:
  ---          if P.F=FConst then P.C(I) is assumed to be zero for I>0

  --------------- standard Scalar procedures

  --- basic
  function Info(Dummy: Taylor1) return Scalar_Info;
  function IsSharp(P: Taylor1) return Boolean;         --- polynomial with sharp coeffs
  function IsZero(P: Taylor1) return Boolean;
  function "="(P1,P2: Taylor1) return Boolean;
  procedure SetZero(P: in out Taylor1);
  procedure Copy(P1: in Taylor1; P2: in out Taylor1);
  procedure Swap(P1,P2: in out Taylor1);

  --- sets
  function Center0(P: Taylor1) return Boolean;         --- symmetric about 0
  function Contains0(P: Taylor1) return Logical;
  function Contains(P1,P2: Taylor1) return Logical;
  procedure BallAt0(R: in Flt; P: in out Taylor1);     --- 0-centered ball of radius |R|
  procedure DiskAt0(F: in Flt; S: in out Taylor1) renames BallAt0;
  function BallAt0(R: Flt) return Taylor1;             --- 0-centered ball of radius |R|
  function DiskAt0(F: Flt) return Taylor1 renames BallAt0;
  procedure ToErr(P: in out Taylor1);                  --- 0-centered set containing P
  procedure ToErr(P1: in Taylor1; P2: in out Taylor1); --- 0-centered set containing P1
  function ToErr(P: Taylor1) return Taylor1;           --- 0-centered set containing P
  procedure Center(P: in out Taylor1);
  procedure Center(P1: in Taylor1; P2: in out Taylor1);
  function Center(P: Taylor1) return Taylor1;
  procedure ModCenter(P: in out Taylor1);
  procedure ModCenter(P1: in Taylor1; P2: in out Taylor1);
  function ModCenter(P: Taylor1) return Taylor1;                      --- zero center
  procedure ErrMult(R: in Radius; P: in out Taylor1);
  procedure Union(P1: in Taylor1; P2: in out Taylor1);
  function Union(P1,P2: Taylor1) return Taylor1;
  procedure Intersection(P1: in Taylor1; P2: in out Taylor1; Empty: out Logical);
  function Interval(C: Flt; R: Radius) return Taylor1;

  --- order
  function Sign(P: Taylor1) return Integer;
  function Compare(P1,P2: Taylor1) return Integer;
  function "<"(P1,P2: Taylor1) return Boolean;
  function "<="(P1,P2: Taylor1) return Boolean;
  function ">="(P1,P2: Taylor1) return Boolean;
  function ">"(P1,P2: Taylor1) return Boolean;
  procedure Min(P1: in Taylor1; P2: in out Taylor1);
  procedure Min(P1,P2: in Taylor1; P3: in out Taylor1);
  function Min(P1,P2: Taylor1) return Taylor1;
  procedure Max(P1: in Taylor1; P2: in out Taylor1);
  procedure Max(P1,P2: in Taylor1; P3: in out Taylor1);
  function Max(P1,P2: Taylor1) return Taylor1;
  function Inf(P: Taylor1; Prec: Radius) return Flt; -- inf of real part for real arguments
  function Inf(P: Taylor1) return Flt;               -- ... with Prec=Psi
  function Sup(P: Taylor1; Prec: Radius) return Flt; -- sup of real part for real arguments
  function Sup(P: Taylor1) return Flt;               -- ... with Prec=Psi

  --- addition and multiplication etc.
  procedure Neg(P: in out Taylor1);
  procedure Neg(P1: in Taylor1; P2: in out Taylor1);
  function "-"(P: Taylor1) return Taylor1;
  procedure Add(I: in Integer; P: in out Taylor1);
  procedure Add(P1: in Taylor1; P2: in out Taylor1);
  procedure Sum(P1,P2: in Taylor1; P3: in out Taylor1);
  function "+"(P1,P2: Taylor1) return Taylor1;
  procedure Sub(P1: in Taylor1; P2: in out Taylor1);
  procedure Diff(P1,P2: in Taylor1; P3: in out Taylor1);
  function "-"(P1,P2: Taylor1) return Taylor1;
  procedure Mult(R: in Flt; P,Tmp: in out Taylor1);
  procedure Mult(R: in Flt; P: in out Taylor1; Tmp: in out Scalar);
  procedure Mult(R: in Flt; P: in out Taylor1);
  procedure Prod(R: in Flt; P1: in Taylor1; P2: in out Taylor1);
  function "*"(R: Flt; P: Taylor1) return Taylor1;
  procedure AddProd(R: in Flt; P1: in Taylor1; P2,Tmp: in out Taylor1);
  procedure AddProd(R: in Flt; P1: in Taylor1; P2: in out Taylor1);
  procedure Mult(Q: in Rational; P: in out Taylor1);
  procedure Mult(P1: in Taylor1; P2,Tmp: in out Taylor1);
  procedure Mult(P1: in Taylor1; P2: in out Taylor1);
  procedure Prod(P1,P2: in Taylor1; P3: in out Taylor1);
  function "*"(P1,P2: Taylor1) return Taylor1;
  procedure AddProd(P1,P2: in Taylor1; P3: in out Taylor1);
  procedure SubProd(P1,P2: in Taylor1; P3: in out Taylor1);
  procedure AddProd(P1,P2: in Taylor1; P3,Dummy: in out Taylor1);
  procedure SubProd(P1,P2: in Taylor1; P3,Dummy: in out Taylor1);
  procedure Div(R: in Flt; P,Tmp: in out Taylor1);
  procedure Div(R: in Flt; P: in out Taylor1; Tmp: in out Scalar);
  procedure Div(R: in Flt; P: in out Taylor1);
  procedure Quot(P1: in Taylor1; R: in Flt; P2: in out Taylor1);
  function "/"(P: Taylor1; R: Flt) return Taylor1;
  procedure Div(P1: in Taylor1; P2: in out Taylor1);
  procedure Quot(P1,P2: in Taylor1; P3: in out Taylor1);
  function "/"(P1,P2: Taylor1) return Taylor1;
  procedure Inv(P: in out Taylor1);
  procedure Inv(P1: in Taylor1; P2: in out Taylor1);
  function Inv(P: Taylor1) return Taylor1;
  function "**"(P: Taylor1; I: Integer) return Taylor1;

  --- functions
  function IsReal(P: Taylor1) return Boolean;
  procedure Adjoint(P: in out Taylor1);
  procedure Adjoint(P1: in Taylor1; P2: in out Taylor1);
  function Adjoint(P: Taylor1) return Taylor1;
  procedure Real_Part(P: in out Taylor1);
  procedure Imag_Part(P: in out Taylor1);
  procedure Eval0(P: in out Taylor1);
  procedure Norm(P: in Taylor1; Q: in out Taylor1);
  function Norm(P: Taylor1) return Taylor1;
  function "abs"(P: Taylor1) return Taylor1 renames Norm;
  function MaxNorm(P: Taylor1) return Radius;
  procedure Sqr(P: in out Taylor1);
  function Sqr(P: Taylor1) return Taylor1;
  procedure Sqr(P1: in Taylor1; P2: in out Taylor1);
  procedure Sqrt(P1: in Taylor1; P2: in out Taylor1);
  function Sqrt(P: Taylor1) return Taylor1;
  procedure Root(K: in Positive; P1: in Taylor1; P2: in out Taylor1);
  function Root(K: Positive; P: Taylor1) return Taylor1;
  procedure Exp(P1: in Taylor1; P2: in out Taylor1);
  function Exp(P: Taylor1) return Taylor1;
  procedure ArcCos(P1: in Taylor1; P2: in out Taylor1);
  function ArcCos(P: Taylor1) return Taylor1;
  procedure ArcSin(P1: in Taylor1; P2: in out Taylor1);
  function ArcSin(P: Taylor1) return Taylor1;
  procedure Cos(P1: in Taylor1; P2: in out Taylor1);
  function Cos(P: Taylor1) return Taylor1;
  procedure Sin(P1: in Taylor1; P2: in out Taylor1);
  function Sin(P: Taylor1) return Taylor1;
  procedure Log(P1: in Taylor1; P2: in out Taylor1);
  function Log(P: Taylor1) return Taylor1;
  procedure Simple_Random(P: in out Taylor1);
  procedure QPower(Q: in Rational; P1: in Taylor1; P2: in out Taylor1);
  procedure Roots2(B,C: in Taylor1; U1,U2,V: in out Taylor1);
  procedure Roots3(B,C,D: in Taylor1; U0,U1,U2,V: in out Taylor1);

  --- conversion and i/o
  function Approx(P: Taylor1) return Flt;
  procedure Assign(I: in Integer; P: in out Taylor1);
  procedure Assign(Q: in Rational; P: in out Taylor1);
  procedure Assign(R: in Flt; P: in out Taylor1);
  function Scal(I: Integer) return Taylor1;
  function Scal(Q: Rational) return Taylor1;
  function Scal(R: Flt) return Taylor1;
  procedure Enclose(R1,R2: in Flt; P: in out Taylor1);
  function Enclose(R1,R2: Flt) return Taylor1;
  procedure Show1(N: in String; P: in Taylor1; Hide0: in Boolean := True);
  procedure Show2(N: in String; P1,P2: in Taylor1; Hide0: in Boolean := True);
  procedure Put(F: in File_Type; P: in Taylor1; Decimal: in Boolean := False);
  procedure Get(F: in File_Type; P: in out Taylor1; Decimal: in Boolean := False);
  procedure Write(FileName: in String; P: in Taylor1; Decimal: in Boolean := False);
  procedure Read(FileName: in String; P: in out Taylor1; Decimal: in Boolean := False);
  function Read(FileName: String; Decimal: Boolean := False) return Taylor1;

  --- misc
  function Get_Precision(P: Taylor1) return Positive;
  procedure Proper_Rounding(Dummy: in Taylor1);
  procedure Free_Cache(Dummy: in Taylor1);

  --------------- other standard procedures
  procedure Show_Param_T1(N: in String := "");
  function Rho(P: Taylor1) return Radius;
  procedure SetRho(R: in Radius; P: in out Taylor1; Check: in Boolean := True);
  procedure AdjustF(P: in out Taylor1);                                -- adjust error index F
  procedure AdjustF(P: in out Taylor1; D: out Natural);                -- adjust and determine EffDeg
  procedure NonConst(P: in out Taylor1);                               -- use before modifying P.C
  procedure LiftErrs(P: in out Taylor1; F: in Power := 0);             -- make general errors start at degree F
  procedure HiErrs(P: in out Taylor1);                                 -- LiftErrs(P,PDeg)
  procedure Assign(S: in Scalar; P: in out Taylor1);
  function Scal(S: Scalar) return Taylor1;
  procedure SetZero(R: in Radius; P: in out Taylor1);
  procedure SetOne(R: in Radius; P: in out Taylor1);
  procedure Affine(C0,C1: in Scalar; P: in out Taylor1);
  procedure Component(I: in Natural; P: in Taylor1; S: in out Scalar);
  procedure Coeff(I: in Natural; P: in Taylor1; S: in out Scalar);
  function Coeff(I: Natural; P: Taylor1) return Scalar;
  procedure AddCoeff(I: in Natural; S: in Scalar; P: in out Taylor1);       -- P[I] := P[I] + S
  procedure Extract_Coeff(I: in Power; P: in out Taylor1; S: in out Scalar);
  procedure Add_Ball(I: in Power; R: in Flt; P: in out Taylor1);

  procedure Add(S: in Scalar; P: in out Taylor1);
  procedure Mult(S: in Scalar; P: in out Taylor1);                      -- P := S*P
  procedure Prod(S: in Scalar; P1: in Taylor1; P2: in out Taylor1);     -- P2 := S*P1;
  function "*"(S: Scalar; P: Taylor1) return Taylor1;                   -- return S*P
  procedure AddProd(S: in Scalar; P1: in Taylor1; P2: in out Taylor1);  -- P2 := P2+S*P1
  procedure SubProd(S: in Scalar; P1: in Taylor1; P2: in out Taylor1);  -- P2 := P2-S*P1
  procedure Div(S: in Scalar; P: in out Taylor1);

  procedure Val(P: in Taylor1; S1: in Scalar; S2: in out Scalar);
  function Val(P: Taylor1; S: Scalar) return Scalar;
  procedure Val(P1: in Taylor1; S: in Scalar; P2: in out Taylor1);
  procedure ValBall(P: in Taylor1; R: in Radius; S: in out Scalar);
  procedure Ran(P: in Taylor1; S: in out Scalar);
  function MaxNorm(R: Radius; P: Taylor1) return Radius;
  function MaxDist(R: Radius; P1,P2: Taylor1) return Radius;
  procedure Norm1(R: in Radius; P: in Taylor1; S: in out Scalar);
  procedure Norm1(P: in Taylor1; S: in out Scalar);
  function Norm1(R: Radius; P: Taylor1) return Scalar;
  function Norm1(P: Taylor1) return Scalar;
  procedure Monom(I: in Power; P: in out Taylor1);
  procedure Compose(C: in SP.Polynom1; P1: in Taylor1; P2: in out Taylor1);
  procedure Reflect(P: in out Taylor1; Offset: in Parity);
  procedure CutParity(P: in out Taylor1; Offset: in Parity);
  function ParityErr(P: Taylor1; Offset: Parity) return Radius;
  procedure Scale(S: in Scalar; P: in out Taylor1);
  procedure Scale(P1: in Taylor1; S: in Scalar; P2: in out Taylor1);
  procedure Translate(P1: in Taylor1; S: in Scalar; P2: in out Taylor1);
  procedure Translate(P1: in Taylor1; S: in Scalar; P2: in out Taylor1; R: in Radius); -- P2 := P1(.+S)
  procedure CosSin(P: in Taylor1; PC,PS: in out Taylor1);

  procedure Der(R: in Radius; P: in out Taylor1);
  procedure ValDer(P: in Taylor1; S: in Scalar; D: in out Scalar);
  procedure DivArg(P: in out Taylor1);
  procedure AntiDer(P: in out Taylor1);
  procedure AntiDer(P1: in Taylor1; P2: in out Taylor1);
  function Integral(P: Taylor1; S1,S2: Scalar) return Scalar;
  procedure Inv_MDer(P: in out Taylor1; CheckConst: in Boolean := True);

  function Pi return Taylor1;

  procedure CosTaylor(S: in Scalar; R: in Radius; P: in out Taylor1); --- Taylor series of X --> cos(S+X)
  procedure CosTaylor(K: in Integer; S: in Scalar; R: in Radius; P: in out Taylor1); --- Taylor series of X --> cos(K*(S+X))
  procedure CosTaylor(K,S: in Scalar; R: in Radius; P: in out Taylor1); --- Taylor series of X --> cos(K*(S+X))
  procedure SinTaylor(S: in Scalar; R: in Radius; P: in out Taylor1); --- Taylor series of X --> sin(S+X)
  procedure SinTaylor(K: in Integer; S: in Scalar; R: in Radius; P: in out Taylor1); --- Taylor series of X --> sin(K*(S+X))
  procedure SinTaylor(K,S: in Scalar; R: in Radius; P: in out Taylor1); --- Taylor series of X --> sin(K*(S+X))

  ---- numerical
  function GuessRho(P: Taylor1) return Radius;
  function GuessRhoInv(P: Taylor1) return Radius;
  procedure NumInterpolate(P0: in Taylor1; P: in out Taylor1; Fac: in Flt);
  procedure NumTrunc(P: in out Taylor1; Deg: in Power := PDeg; Cut: in Radius := Zero);
  procedure NumIProd(P1,P2: in Taylor1; S: in out Scalar);

  package TAlloc is new Recycling (Item => Taylor1);
  subtype LT_Pointer is TAlloc.LD_Pointer;
  TAC: TAlloc.Controlled_Protected_Data;
  TPool: TAlloc.Protected_Data renames TAC.P;

private

  function NotConst(C: Poly1) return Boolean;
  procedure AddProd(S: in Scalar; C1: in Poly1; C2: in out Poly1; Tmp: in out Scalar);
  procedure SubProd(S: in Scalar; C1: in Poly1; C2: in out Poly1; Tmp: in out Scalar);
  function QuasiDeg(P: Taylor1) return Natural;
  function GuessMaxAbs(P: Taylor1; D: Natural) return Radius;
  function Sharp_Sup(P: Taylor1; Prec: Radius) return Flt;
  function Simple(P: Taylor1) return Boolean;
  function TrueF(P: Taylor1) return Natural;
  procedure Add_Coeff(I: in Power; S: in Scalar; P: in out Taylor1);
  procedure Sub_Coeff(I: in Power; S: in Scalar; P: in out Taylor1);
  function NotSharp(C: Poly1; F,L: Integer) return Boolean;

  function Prod_Err(C1,C2: Poly1; D1,D2: Natural; R: Radius) return Radius;
  procedure Add_Prod(P,Q: in Taylor1; H: in out Taylor1);
  procedure Num_ValDer(P: in Poly1; S: in Scalar; C: in out Scalar);
  procedure IPower(I: in Integer; P1: in Taylor1; P2: in out Taylor1);
  procedure Inv_Split(P1: in Taylor1; P2: in out Taylor1; S: in out Scalar);
  procedure QPower_Split(Q: in Rational; P1: in Taylor1; P2: in out Taylor1; S: in out Scalar);
  procedure Exp_Split(P1: in Taylor1; P2: in out Taylor1; S: in out Scalar);
  procedure CosSin_Split(P1: in Taylor1; P2: in out Taylor1; C,S: in out Scalar);
  procedure Log_Split(P1: in Taylor1; P2: in out Taylor1; S: in out Scalar);
  procedure ArcSin_Part(P1: in Taylor1; P2: in out Taylor1);
  function ExpTail(K: Natural; R: Radius) return Radius;

  PDeg1:   constant Integer     := PDeg+1;
  SZero:   constant Scalar      := Scal(0);
  SOne:    constant Scalar      := Scal(1);
  SInfo:   constant Scalar_Info := Info(SZero);
  STrunc:  constant Boolean     := SInfo.IsNumeric;
  Not_STrunc: constant Boolean  := not STrunc;
  Std:     constant Boolean     := SInfo.IsStandard;
  TInfo:   constant Scalar_Info := (RepType    => SInfo.RepType,
                                    IsStandard => SInfo.IsStandard,
                                    IsNumeric  => SInfo.IsNumeric,
                                    IsComplex  => SInfo.IsComplex,
                                    NReps      => SInfo.NReps*PDeg1,
                                    RMode      => SInfo.RMode);

  pragma Inline_Always (IsZero,IsSharp,Simple);

  -------------- FBall

  type FBall is
    record
      C: Flt;    -- center
      R: Radius; -- radius
    end record;

  type FPoly1 is array(Power) of FBall;

  function MaxNorm(B: in FBall) return Radius;
  procedure Enclose(R1,R2: in Flt; B: out FBall);
  procedure AddProd(B1,B2: in FBall; B3: in out FBall);

  package BAlloc is new Recycling (Item => FPoly1);
  subtype LB_Pointer is BAlloc.LD_Pointer;
  BAC: BAlloc.Controlled_Protected_Data;
  BPool: BAlloc.Protected_Data renames BAC.P;

end Taylors1;
