with Strings, Globals, Ints, IPowers, Fun_Series, Newton, Roots;
use Strings, Globals, Ints;

pragma Elaborate_All (Strings,Globals,Ints,IPowers,Fun_Series,Newton, Roots);

package body Taylors2 is

  procedure SubProd(S: in Scalar; C1: in Poly2; C2: in out Poly2) is
  begin
    if not IsZero(S) then
      declare
        Tmp: Scalar;
      begin
        for I in Power loop
          for J in 0 .. PDeg-I loop SubProd(S,C1(I,J),C2(I,J),Tmp); end loop;
        end loop;
      end;
    end if;
  end SubProd;

  procedure Add(C: in Poly2C; Deg: in Natural; P: in out Taylor2) is
    PC: Poly2 renames P.C;
    PR: Radius renames P.R;
    N,J: Integer;
    Fac: Flt;
    Tmp: Scalar renames PC(PDeg,PDeg);
    E: Scalar;
  begin
    Add(C,PC); -- coeffs up do degree PDeg
    for D in reverse PDeg1 .. Deg loop
      N := D-PDeg;
      Fac := PR**N;
      for I in 0 .. D loop
        J := D-I;
        ToErr(C(I,J),E);
        if I>J then
          AddProd(Fac,E,PC(I-N,J),Tmp);
        elsif J>I then
          AddProd(Fac,E,PC(I,J-N),Tmp);
        else
          AddProd(Half*Fac,E,PC(I-N,J),Tmp);
          AddProd(Half*Fac,E,PC(I,J-N),Tmp);
        end if;
      end loop;
    end loop;
  end Add;

  --- basic

  function Info(Dummy: Taylor2) return Scalar_Info is
    pragma Unreferenced(Dummy);
  begin
    return TInfo;
  end Info;

  function IsSharp(P: Taylor2) return Boolean is
  begin
    return STrunc or else IsSharp(P.C);
  end IsSharp;

  function IsZero(P: Taylor2) return Boolean is
  begin
    if P.F=FConst then
      return IsZero(P.C(0,0));
    else
      return IsZero(P.C);
    end if;
  end IsZero;

  function "="(P1,P2: Taylor2) return Boolean is
    C1: Poly2 renames P1.C;
    C2: Poly2 renames P2.C;
  begin
    if (C1 /= C2) then return False; end if;
    if IsSharp(P1) then return True; end if;
    if (P1.R /= P2.R) then return False; end if;
    return TrueF(P1)=TrueF(P2);
  end "=";

  procedure SetZero(P: in out Taylor2) is
  begin
    SetZero(P.C);
    P.F := FConst;
    P.R := Rad;
  end SetZero;

  procedure Copy(P1: in Taylor2; P2: in out Taylor2) is
  begin
    Copy(P1.C,P2.C);
    P2.F := P1.F;
    P2.R := P1.R;
  end Copy;

  procedure Swap(P1,P2: in out Taylor2) is
    F1: constant Integer := P1.F;
    R1: constant Flt := P1.R;
  begin
    Swap(P1.C,P2.C);
    P1.F := P2.F; P2.F := F1;
    P1.R := P2.R; P2.R := R1;
  end Swap;

  --- sets

  function Center0(P: Taylor2) return Boolean is
  begin
    return Center0(P.C);
  end Center0;

  function Contains0(P: Taylor2) return Logical is
    --- very crude
    L: Logical := Contains0(P.C(0,0));
  begin
    if P.F=FConst then return L; end if;
    if L=False then return False; end if;
    L := Contains0(P.C);
    if L=True or else STrunc then return L; end if;
    return Uncertain;
  end Contains0;

  function Contains(P1,P2: Taylor2) return Logical is
    --- very crude
    C1: Poly2 renames P1.C;
    C2: Poly2 renames P2.C;
    L: Logical := Contains(C1(0,0),C2(0,0));
  begin
    if L=False then return False; end if;
    L := Contains(C1,C2);
    if STrunc then return L; end if;
    if L=True then
      if P2.F>PDeg then return True; end if;
      if (P1.R>P2.R) then return False; end if;
      if P1.F=PDeg and then P2.F=PDeg then return True; end if;
      declare
        D1: constant Integer := TrueF(P1);
        D2: constant Integer := TrueF(P2);
      begin
        if D1=PDeg and then D2=PDeg then return True; end if;
        if D1>D2 then return False; end if;
      end;
    end if;
    return Uncertain;
  end Contains;

  procedure BallAt0(R: in Flt; P: in out Taylor2) is
  begin
    SetZero(P.C);
    if STrunc then
      P.F := FConst;
    else
      BallAt0(R,P.C(0,0));
      P.F := 0;
    end if;
    P.R := Rad;
  end BallAt0;

  function BallAt0(R: Flt) return Taylor2 is
    P: Taylor2;
  begin
    BallAt0(R,P);
    return P;
  end BallAt0;

  procedure DiskAt0(R: in Flt; P: in out Taylor2) is
  begin
    SetZero(P.C);
    P.F := FConst;
    if not STrunc then
      BallAt0(R,P.C(0,0));
    end if;
    P.R := Rad;
  end DiskAt0;

  function DiskAt0(R: Flt) return Taylor2 is
    P: Taylor2;
  begin
    DiskAt0(R,P);
    return P;
  end DiskAt0;

  procedure ToErr(P: in out Taylor2) is
  begin
    P.F := TrueF(P);
    ToErr(P.C);
  end ToErr;

  procedure ToErr(P1: in Taylor2; P2: in out Taylor2) is
  begin
    ToErr(P1.C,P2.C);
    P2.F := TrueF(P1);
    P2.R := P1.R;
  end ToErr;

  function ToErr(P: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    ToErr(P,Q);
    return Q;
  end ToErr;

  procedure Center(P: in out Taylor2) is
  begin
    if not STrunc then
      Center(P.C);
      if QuasiDeg0(P.C) then
        P.F := FConst;
      else
        P.F := PDeg1;
      end if;
    end if;
  end Center;

  procedure Center(P1: in Taylor2; P2: in out Taylor2) is
  begin
    Center(P1.C,P2.C);
    if QuasiDeg0(P2.C) then
      P2.F := FConst;
    else
      P2.F := PDeg1;
    end if;
    P2.R := Rad;
  end Center;

  function Center(P: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    Center(P,Q);
    return Q;
  end Center;

  procedure ModCenter(P: in out Taylor2) is
  begin
    if STrunc then
      SetZero(P);
    else
      ModCenter(P.C);
      P.F := TrueF(P);
    end if;
  end ModCenter;

  procedure ModCenter(P1: in Taylor2; P2: in out Taylor2) is
  begin
    if STrunc then
      SetZero(P2);
    else
      ModCenter(P1.C,P2.C);
      P2.F := TrueF(P2);
      P2.R := P1.R;
    end if;
  end ModCenter;

  function ModCenter(P: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    ModCenter(P,Q);
    return Q;
  end ModCenter;

  procedure ErrMult(R: in Radius; P: in out Taylor2) is
  begin
    if not STrunc then ErrMult(R,P.C); end if;
  end ErrMult;

  procedure Union(P1: in Taylor2; P2: in out Taylor2) is
  begin
    Union(P1.C,P2.C);
    if P2.F<P1.F then P2.F := P1.F; end if;
    if P2.R>P1.R then P2.R := P1.R; end if;
  end Union;

  function Union(P1,P2: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    Copy(P2,Q);
    Union(P1,Q);
    return Q;
  end Union;

  procedure Intersection(P1: in Taylor2; P2: in out Taylor2; Empty: out Logical) is
    --- far from optimal
    C1: Poly2 renames P1.C;
    C2: Poly2 renames P2.C;
    F1: constant Integer := TrueF(P1);
    F2: Natural renames P2.F;
    E: Logical;
    Tmp: LT_Pointer;
  begin
    F2 := TrueF(P2);
    if F1>F2 then
      Pool.Allocate(Tmp);
      Copy(P2,Tmp.all.Data);
      Copy(P1,P2);
      Intersection(Tmp.all.Data,P2,Empty);
      Pool.Recycle(Tmp);
    else
      if P2.R<P1.R then P2.R := P1.R; end if;
      for D in 0 .. IMin(F1,PDeg) loop
        for I in 0 .. D loop
          Intersection(C1(I,D-I),C2(I,D-I),E);
          Empty := Empty or E;
        end loop;
      end loop;
      if F1<PDeg then Empty := Empty or Uncertain; end if;
    end if;
  end Intersection;

  --- order

  function VarNorm(P: Taylor2) return Radius is
    --- assumes Flt operations are rounded up
    C: Poly2 renames P.C;
    R: constant Radius := P.R;
    N: Flt := Zero;
  begin
    for D in reverse 1 .. QuasiDeg(C) loop
      for I in 0 .. D loop N := N+MaxNorm(C(I,D-I)); end loop;
      N := R*N;
    end loop;
    return N;
  end VarNorm;

  function Ran(P: Taylor2) return Scalar is
    --- primitive bound on the range
  begin
    if P.F=FConst then return P.C(0,0); end if;
    declare
      S: Scalar;
    begin
      BallAt0(VarNorm(P),S);
      Add(P.C(0,0),S);
      return S;
    end;
  end Ran;

  function Sign(P: Taylor2) return Integer is
  begin
    if not IsReal(P) then
      raise Undefined with Show0("Taylor2.Sign error: argument non-real");
    end if;
    return Sign(Ran(P));
  end Sign;

  function Simple(P: Taylor2) return Boolean is
    --- only P.C(0,0) nonzero
  begin
    return (P.F=FConst) or else QuasiDeg0(P.C);
  end Simple;

  function Compare(P1,P2: Taylor2) return Integer is
  begin
    if Simple(P1) and then Simple(P2) then
      return Compare(P1.C(0,0),P2.C(0,0));
    else
      declare
        S: Integer;
        Tmp: LT_Pointer;
      begin
        Pool.Allocate(Tmp);
        Diff(P1,P2,Tmp.all.Data);
        S := Sign(Tmp.all.Data);
        Pool.Recycle(Tmp);
        return S;
      end;
    end if;
  end Compare;

  function "<"(P1,P2: Taylor2) return Boolean is
  begin
    return (Compare(P1,P2) < 0);
  end "<";

  function "<="(P1,P2: Taylor2) return Boolean is
  begin
    return (Compare(P1,P2) <= 0);
  end "<=";

  function ">="(P1,P2: Taylor2) return Boolean is
  begin
    return (Compare(P1,P2) >= 0);
  end ">=";

  function ">"(P1,P2: Taylor2) return Boolean is
  begin
    return (Compare(P1,P2) > 0);
  end ">";

  procedure Min(P1: in Taylor2; P2: in out Taylor2) is
  begin
    if P1<P2 then Copy(P1,P2); end if;
  end Min;

  procedure Min(P1,P2: in Taylor2; P3: in out Taylor2) is
  begin
    if P1<P2 then Copy(P1,P3); else Copy(P2,P3); end if;
  end Min;

  function Min(P1,P2: Taylor2) return Taylor2 is
  begin
    if P1<P2 then return P1; else return P2; end if;
  end Min;

  procedure Max(P1: in Taylor2; P2: in out Taylor2) is
  begin
    if P1>P2 then Copy(P1,P2); end if;
  end Max;

  procedure Max(P1,P2: in Taylor2; P3: in out Taylor2) is
  begin
    if P1>P2 then Copy(P1,P3); else Copy(P2,P3); end if;
  end Max;

  function Max(P1,P2: Taylor2) return Taylor2 is
  begin
    if P1>P2 then return P1; else return P2; end if;
  end Max;

  function Sup(P: Taylor2) return Flt is
  begin
    if not IsReal(P) then
      raise Undefined with Show0("Taylor2.Sup error: argument non-real");
    end if;
    return Sup(Ran(P));
  end Sup;

  function Inf(P: Taylor2) return Flt is
  begin
    if not IsReal(P) then
      raise Undefined with Show0("Taylor2.Inf error: argument non-real");
    end if;
    return Inf(Ran(P));
  end Inf;

  --- addition and multiplication etc.

  procedure Neg(P: in out Taylor2) is
  begin
    Neg(P.C);
  end Neg;

  procedure Neg(P1: in Taylor2; P2: in out Taylor2) is
  begin
    Neg(P1.C,P2.C);
    P2.F := P1.F;
    P2.R := P1.R;
  end Neg;

  function "-"(P: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    Neg(P,Q);
    return Q;
  end "-";

  procedure Add(I: in Integer; P: in out Taylor2) is
  begin
    Add(I,P.C(0,0));
  end Add;

  procedure Add(P1: in Taylor2; P2: in out Taylor2) is
  begin
    Add(P1.C,P2.C);
    if P2.F>P1.F then P2.F := P1.F; end if;
    if P2.R>P1.R then P2.R := P1.R; end if;
  end Add;

  procedure Sum(P1,P2: in Taylor2; P3: in out Taylor2) is
  begin
    Sum(P1.C,P2.C,P3.C);
    P3.F := IMin(P1.F,P2.F);
    P3.R := RMin(P1.R,P2.R);
  end Sum;

  function "+"(P1,P2: Taylor2) return Taylor2 is
    P3: Taylor2;
  begin
    Sum(P1,P2,P3);
    return P3;
  end "+";

  procedure Sub(P1: in Taylor2; P2: in out Taylor2) is
  begin
    Sub(P1.C,P2.C);
    if P2.F>P1.F then P2.F := P1.F; end if;
    if P2.R>P1.R then P2.R := P1.R; end if;
  end Sub;

  procedure Diff(P1,P2: in Taylor2; P3: in out Taylor2) is
  begin
    Diff(P1.C,P2.C,P3.C);
    P3.F := IMin(P1.F,P2.F);
    P3.R := RMin(P1.R,P2.R);
  end Diff;

  function "-"(P1,P2: Taylor2) return Taylor2 is
    P3: Taylor2;
  begin
    Diff(P1,P2,P3);
    return P3;
  end "-";

  procedure Mult(R: in Flt; P,Tmp: in out Taylor2) is
    STmp: Scalar renames Tmp.C(0,0);
  begin
    if R=Zero then
      SetZero(P);
    else
      Mult(R,P.C,STmp);
    end if;
  end Mult;

  procedure Mult(R: in Flt; P: in out Taylor2) is
  begin
    if R=Zero then
      SetZero(P);
    else
      declare
        STmp: Scalar;
      begin
        Mult(R,P.C,STmp);
      end;
    end if;
  end Mult;

  procedure Prod(R: in Flt; P1: in Taylor2; P2: in out Taylor2) is
  begin
    if R=Zero then
      SetZero(P2);
    else
      Prod(R,P1.C,P2.C);
      P2.F := P1.F;
      P2.R := P1.R;
    end if;
  end Prod;

  function "*"(R: Flt; P: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    Prod(R,P,Q);
    return Q;
  end "*";

  procedure AddProd(R: in Flt; P1: in Taylor2; P2,Tmp: in out Taylor2) is
    STmp: Scalar renames Tmp.C(0,0);
  begin
    if R /= Zero then
      AddProd(R,P1.C,P2.C,STmp);
      if P2.F>P1.F then P2.F := P1.F; end if;
      if P2.R>P1.R then P2.R := P1.R; end if;
    end if;
  end AddProd;

  procedure AddProd(R: in Flt; P1: in Taylor2; P2: in out Taylor2) is
  begin
    if R /= Zero then
      declare
        STmp: Scalar;
      begin
        AddProd(R,P1.C,P2.C,STmp);
      end;
      if P2.F>P1.F then P2.F := P1.F; end if;
      if P2.R>P1.R then P2.R := P1.R; end if;
    end if;
  end AddProd;

  procedure Mult(Q: in Rational; P: in out Taylor2) is
  begin
    if INum(Q)=0 then
      SetZero(P);
    else
      declare
        Tmp: Scalar;
      begin
        Mult(Flt(INum(Q)),P.C,Tmp);
        Div(Flt(IDen(Q)),P.C,Tmp);
      end;
    end if;
  end Mult;

  procedure Mult(P1: in Taylor2; P2,Tmp: in out Taylor2) is
  begin
    if P1.F=FConst then
      Mult(P1.C(0,0),P2);
    elsif P2.F=FConst then
      declare
        STmp: Scalar renames Tmp.C(0,0);
      begin
        Copy(P2.C(0,0),STmp);
        Prod(STmp,P1,P2);
      end;
    else
      Copy(P2,Tmp);
      SetZero(P2);
      AddProd(P1,Tmp,P2);
    end if;
  end Mult;

  procedure Mult(P1: in Taylor2; P2: in out Taylor2) is
    PTmp: LT_Pointer;
  begin
    if P1.F=FConst then
      Mult(P1.C(0,0),P2);
    elsif P2.F=FConst then
      declare
        STmp: Scalar;
      begin
        Copy(P2.C(0,0),STmp);
        Prod(STmp,P1,P2);
      end;
    else
      Pool.Allocate(PTmp);
      Copy(P2,PTmp.all.Data);
      SetZero(P2);
      AddProd(P1,PTmp.all.Data,P2);
      Pool.Recycle(PTmp);
    end if;
  end Mult;

  procedure Prod(P1,P2: in Taylor2; P3: in out Taylor2) is
  begin
    if P1.F=FConst then
      Prod(P1.C(0,0),P2,P3);
    elsif P2.F=FConst then
      Prod(P2.C(0,0),P1,P3);
    else
      SetZero(P3);
      AddProd(P1,P2,P3);
    end if;
  end Prod;

  function "*"(P1,P2: Taylor2) return Taylor2 is
    P3: Taylor2;
  begin
    Prod(P1,P2,P3);
    return P3;
  end "*";

  procedure AddProd(P1,P2: in Taylor2; P3,Dummy: in out Taylor2) is
    pragma Unreferenced (Dummy);
  begin
    AddProd(P1,P2,P3);
  end AddProd;

  procedure AddProd(P1,P2: in Taylor2; P3: in out Taylor2) is
    C1: Poly2 renames P1.C;
    C2: Poly2 renames P2.C;
  begin
    if P1.F=FConst then
      AddProd(C1(0,0),P2,P3);
      P3.F := IMin(P2.F,P3.F);
      P3.R := RMin(P2.R,P3.R);
    elsif P2.F=FConst then
      AddProd(C2(0,0),P1,P3);
      P3.F := IMin(P1.F,P3.F);
      P3.R := RMin(P1.R,P3.R);
    else
      declare
        D1: constant Integer := QuasiDeg(C1);
        D2: constant Integer := QuasiDeg(C2);
        D3: Integer := D1+D2;
        C3: Poly2 renames P3.C;
        Tmp: LT_PointerC;
      begin
        P3.R := RMin(P1.R,P2.R,P3.R);
        P3.F := IMin(P1.F,P2.F,P3.F);
        if D1=0 then AddProd(C1(0,0),P2,P3); return; end if;
        if D2=0 then AddProd(C2(0,0),P1,P3); return; end if;
        if STrunc or else D3 <= PDeg then
          D3 := IMin(D3,PDeg);
          Add_Prod(C1,C2,D1,D2,D3,C3);
        else
          PoolC.Allocate(Tmp);
          SetZero(Tmp.all.Data,D3);
          Add_Prod(C1,C2,D1,D2,D3,Tmp.all.Data);
          Add(Tmp.all.Data,D3,P3);
          PoolC.Recycle(Tmp);
        end if;
      end;
    end if;
  end AddProd;

  procedure SubProd(P1,P2: in Taylor2; P3,Dummy: in out Taylor2) is
    pragma Unreferenced (Dummy);
  begin
    Neg(P3);
    AddProd(P1,P2,P3);
    Neg(P3);
  end SubProd;

  procedure Div(R: in Flt; P,Tmp: in out Taylor2) is
    STmp: Scalar renames Tmp.C(0,0);
  begin
    Div(R,P.C,STmp);
  end Div;

  procedure Div(R: in Flt; P: in out Taylor2) is
    STmp: Scalar;
  begin
    Div(R,P.C,STmp);
  end Div;

  procedure Quot(P1: in Taylor2; R: in Flt; P2: in out Taylor2) is
  begin
    Quot(P1.C,R,P2.C);
    P2.F := P1.F;
    P2.R := P1.R;
  end Quot;

  function "/"(P: Taylor2; R: Flt) return Taylor2 is
    Q: Taylor2;
  begin
    Quot(P,R,Q);
    return Q;
  end "/";

  procedure Div(P1: in Taylor2; P2: in out Taylor2) is
    Tmp: LT_Pointer;
  begin
    Pool.Allocate(Tmp);
    Inv(P1,Tmp.all.Data);
    Mult(Tmp.all.Data,P2);
    Pool.Recycle(Tmp);
  end Div;

  procedure Quot(P1,P2: in Taylor2; P3: in out Taylor2) is
    Tmp: LT_Pointer;
  begin
    Pool.Allocate(Tmp);
    Inv(P2,Tmp.all.Data);
    SetZero(P3);
    AddProd(P1,Tmp.all.Data,P3);
    Pool.Recycle(Tmp);
  end Quot;

  function "/"(P1,P2: Taylor2) return Taylor2 is
    P3: Taylor2;
  begin
    Quot(P1,P2,P3);
    return P3;
  end "/";

  procedure Inv(P: in out Taylor2) is
    Tmp: LT_Pointer;
  begin
    Pool.Allocate(Tmp);
    Inv(P,Tmp.all.Data);
    Copy(Tmp.all.Data,P);
    Pool.Recycle(Tmp);
  end Inv;

  function Inv(P: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    Inv(P,Q);
    return Q;
  end Inv;

  --- functions

  function IsReal(P: Taylor2) return Boolean is
  begin
    if TInfo.IsComplex then
      if P.F=FConst then
        return IsReal(P.C(0,0));
      else
        return IsReal(P.C);
      end if;
    else
      return True;
    end if;
  end IsReal;

  procedure Adjoint(P: in out Taylor2) is
  begin
    if TInfo.IsComplex then Adjoint(P.C); end if;
  end Adjoint;

  procedure Adjoint(P1: in Taylor2; P2: in out Taylor2) is
  begin
    Adjoint(P1.C,P2.C);
    P2.F := P1.F;
    P2.R := P1.R;
  end Adjoint;

  function Adjoint(P: Taylor2) return Taylor2 is
  begin
    if TInfo.IsComplex then
      declare
        Q: Taylor2;
      begin
        Adjoint(P,Q);
        return Q;
      end;
    end if;
    return P;
  end Adjoint;

  procedure Real_Part(P: in out Taylor2) is
  begin
    if TInfo.IsComplex then
      Real_Part(P.C);
      P.F := TrueF(P);
    end if;
  end Real_Part;

  procedure Imag_Part(P: in out Taylor2) is
  begin
    if TInfo.IsComplex then
      Imag_Part(P.C);
      P.F := TrueF(P);
    else
      SetZero(P);
    end if;
  end Imag_Part;

  procedure Eval0(P: in out Taylor2) is
    C: Poly2 renames P.C;
  begin
    Eval0(C(0,0));
    for D in 1 .. PDeg loop
      for I in 0 .. D loop SetZero(C(I,D-I)); end loop;
    end loop;
    P.F := FConst;
  end Eval0;

  procedure Norm(P: in Taylor2; Q: in out Taylor2) is
  begin
    SetZero(Q);
    Norm1(P,Q.C(0,0));
  end Norm;

  function Norm(P: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    SetZero(Q);
    Norm1(P,Q.C(0,0));
    return Q;
  end Norm;

  function MaxNorm(P: Taylor2) return Radius is
  begin
    if P.F=FConst then
      return MaxNorm(P.C(0,0));
    else
      return MaxNorm(P.C(0,0))+VarNorm(P);
    end if;
  end MaxNorm;

  procedure Sqr(P: in out Taylor2) is
    Tmp: LT_Pointer;
  begin
    if P.F=FConst then
      Sqr(P.C(0,0));
    else
      Pool.Allocate(Tmp);
      Copy(P,Tmp.all.Data);
      SetZero(P);
      AddProd(Tmp.all.Data,Tmp.all.Data,P);
      Pool.Recycle(Tmp);
    end if;
  end Sqr;

  function Sqr(P: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    SetZero(Q);
    AddProd(P,P,Q);
    return Q;
  end Sqr;

  procedure Sqr(P1: in Taylor2; P2: in out Taylor2) is
  begin
    SetZero(P1.R,P2);
    if P1.F=FConst then
      Sqr(P1.C(0,0),P2.C(0,0));
    else
      AddProd(P1,P1,P2);
    end if;
  end Sqr;

  function Sqrt(P: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    Sqrt(P,Q);
    return Q;
  end Sqrt;

  function Root(K: Positive; P: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    Root(K,P,Q);
    return Q;
  end Root;

  function Exp(P: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    Exp(P,Q);
    return Q;
  end Exp;

  function ArcCos(P: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    ArcCos(P,Q);
    return Q;
  end ArcCos;

  function ArcSin(P: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    ArcSin(P,Q);
    return Q;
  end ArcSin;

  procedure Cos(P1: in Taylor2; P2: in out Taylor2) is
    Tmp: LT_Pointer;
  begin
    if P1.F=FConst then
      SetZero(P2);
      Cos(P1.C(0,0),P2.C(0,0));
    else
      Pool.Allocate(Tmp);
      CosSin(P1,P2,Tmp.all.Data);
      Pool.Recycle(Tmp);
    end if;
  end Cos;

  function Cos(P: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    Cos(P,Q);
    return Q;
  end Cos;

  procedure Sin(P1: in Taylor2; P2: in out Taylor2) is
    Tmp: LT_Pointer;
  begin
    if P1.F=FConst then
      SetZero(P2);
      Sin(P1.C(0,0),P2.C(0,0));
    else
      Pool.Allocate(Tmp);
      CosSin(P1,Tmp.all.Data,P2);
      Pool.Recycle(Tmp);
    end if;
  end Sin;

  function Sin(P: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    Sin(P,Q);
    return Q;
  end Sin;

  function Log(P: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    Log(P,Q);
    return Q;
  end Log;

  procedure Simple_Random(P: in out Taylor2) is
    R: constant Radius := P.R*Flt(9)/Flt(8); -- ad hoc
    C: Poly2 renames P.C;
    PowR,Fac: Flt := One;
  begin
    Simple_Random(C(0,0));
    for D in 1 .. PDeg loop
      PowR := PowR/R;
      Fac := PowR/Flt(D+1);
      for I in 0 .. D loop
        Simple_Random(C(I,D-I));
        Mult(Fac,C(I,D-I));
      end loop;
    end loop;
    P.F := PDeg1;
  end Simple_Random;

  --- conversion and i/o

  function Approx(P: Taylor2) return Flt is
  begin
    return Approx(P.C(0,0));
  end Approx;

  procedure Assign(I: in Integer; P: in out Taylor2) is
  begin
    SetZero(P);
    Assign(I,P.C(0,0));
  end Assign;

  procedure Assign(Q: in Rational; P: in out Taylor2) is
  begin
    SetZero(P);
    Assign(Q,P.C(0,0));
  end Assign;

  procedure Assign(R: in Flt; P: in out Taylor2) is
  begin
    SetZero(P);
    Assign(R,P.C(0,0));
  end Assign;

  function Scal(I: Integer) return Taylor2 is
    P: Taylor2;
  begin
    SetZero(P);
    Assign(I,P.C(0,0));
    return P;
  end Scal;

  function Scal(Q: Rational) return Taylor2 is
    P: Taylor2;
  begin
    SetZero(P);
    Assign(Q,P.C(0,0));
    return P;
  end Scal;

  function Scal(R: Flt) return Taylor2 is
    P: Taylor2;
  begin
    SetZero(P);
    Assign(R,P.C(0,0));
    return P;
  end Scal;

  procedure Enclose(R1,R2: in Flt; P: in out Taylor2) is
  begin
    SetZero(P);
    Enclose(R1,R2,P.C(0,0));
  end Enclose;

  function Enclose(R1,R2: Flt) return Taylor2 is
    P: Taylor2;
  begin
    SetZero(P);
    Enclose(R1,R2,P.C(0,0));
    return P;
  end Enclose;

  procedure Show1(N: in String; P: in Taylor2; Hide0: in Boolean := True) is
  begin
    if not IsZero(P) then
      Show1(N & "F ",P.F);
      Show1(N & "R ",P.R);
      Show1(N & "C ",P.C);
    elsif not Hide0 then
      Show0(N & "0");
    end if;
  end Show1;

  procedure Show2(N: in String; P1,P2: in Taylor2; Hide0: in Boolean := True) is
  begin
    if not (IsZero(P1) and IsZero(P2)) then
      Show2(N & "F ",P1.F,P2.F);
      Show2(N & "R ",P1.R,P2.R);
      Show2(N & "C ",P1.C,P2.C);
    elsif not Hide0 then
      Show0(N & "0 0");
    end if;
  end Show2;

  procedure Put(F: in File_Type; P: in Taylor2; Decimal: in Boolean := False) is
    C: Poly2 renames P.C;
    Deg: constant Integer := QuasiDeg(C);
    PF:  constant Integer := IMin(TrueF(P),Deg+1);
  begin
    Put(F,PF);
    Put(F,P.R,Decimal);
    Put(F,Deg);
    for D in 0 .. Deg loop
      for I in 0 .. D loop
        Put(F,P.C(I,D-I),Decimal);
      end loop;
    end loop;
  end Put;

  procedure Get(F: in File_Type; P: in out Taylor2; Decimal: in Boolean := False) is
    Deg: Integer := 0;
    S: Scalar;
  begin
    SetZero(P);
    Get(F,P.F);
    Get(F,P.R,Decimal);
    Get(F,Deg);
    for D in 0 .. Deg loop
      for I in 0 .. D loop
        Get(F,S,Decimal);
        AddCoeff(I,D-I,S,P);
      end loop;
    end loop;
    P.F := TrueF(P);
  end Get;

  procedure Write(FileName: in String; P: in Taylor2; Decimal: in Boolean := False) is
    F: File_Type;
  begin
    if Verbosity>0 then Show0("Writing " & FileName); end if;
    Create(F,Out_File,FileName);
    Put(F,P,Decimal);
    Close(F);
  end Write;

  procedure Read(FileName: in String; P: in out Taylor2; Decimal: in Boolean := False) is
    F: File_Type;
  begin
    if Verbosity>0 then Show0("Reading " & FileName); end if;
    Open(F,In_File,FileName);
    Get(F,P,Decimal);
    Close(F);
  end Read;

  function Read(FileName: String; Decimal: Boolean := False) return Taylor2 is
    P: Taylor2;
  begin
    Read(FileName,P,Decimal);
    return P;
  end Read;

  --- misc

  function Get_Precision(P: Taylor2) return Positive is
  begin
    return Get_Precision(P.C(0,0));
  end Get_Precision;

  procedure Free_Cache(Dummy: in Taylor2) is
    pragma Unreferenced (Dummy);
  begin
    Pool.Free_All;
  end Free_Cache;

  procedure Proper_Rounding(Dummy: in Taylor2) is
    pragma Unreferenced(Dummy);
  begin
    Proper_Rounding;
  end Proper_Rounding;

  --------------- other standard procedures

  function Rho(P: Taylor2) return Radius is
  begin
    return P.R;
  end Rho;

  procedure SetRho(R: in Radius; P: in out Taylor2; Check: in Boolean := True) is
  begin
    P.F := TrueF(P);
    if (P.F <= PDeg) and then Check and then (R>P.R) then
      raise Domain_Error with Show0("Taylors2.SetRho: cannot increase domain");
    end if;
    P.R := R;
  end SetRho;

  function TrueF(P: Taylor2) return Natural is
    --- the value FConst means that P is constant
    --- the value PDeg1 means that P is a polynomial
    C: Poly2 renames P.C;
    Deg: constant Integer := QuasiDeg(C);
  begin
    if STrunc then
      if Deg=0 then return FConst; end if;
      return PDeg1;
    end if;
    if Deg=0 then
      if P.F>0 or else IsSharp(C(0,0)) then return FConst; end if;
      return 0;
    end if;
    for D in P.F .. Deg loop
      for I in 0 .. D loop
        if not IsSharp(C(I,D-I)) then return D; end if;
      end loop;
    end loop;
    return PDeg1;
  end TrueF;

  procedure AdjustF(P: in out Taylor2) is
  begin
    P.F := TrueF(P);
  end AdjustF;

  procedure NonConst(P: in out Taylor2) is
  begin
    if P.F=FConst then P.F := PDeg1; end if;
  end NonConst;

  procedure LiftErrs(P: in out Taylor2; F: in Power := 0) is
  begin
    P.F := TrueF(P);
    if F>P.F then
      declare
        R: Radius renames P.R;
        C: Poly2 renames P.C;
        N: Integer;
        E,Fac: Scalar;
      begin
        if F=1 then
          ModCenter(C(0,0),E);
          Div(R,E);
          Add(E,C(1,0));
          Add(E,C(0,1));
        else
          BallAt0(One,Fac);
          for D in reverse P.F .. F-1 loop -- distribute: degree >= D --> degree >= F
            Div(R,Fac);
            for M in 0 .. D loop
              N := D-M;
              Prod(Fac,C(M,N),E);          -- error from C(M,N) with M+N=D
              for I in M .. F-N loop
                Add(E,C(I,F-I));           -- add to C(I,J) with I+J=F
              end loop;
            end loop;
          end loop;
        end if;
      end;
      P.F := F;
    end if;
  end LiftErrs;

  procedure HiErrs(P: in out Taylor2) is
  begin
    LiftErrs(P,PDeg);
  end HiErrs;

  procedure Assign(S: in Scalar; P: in out Taylor2) is
  begin
    SetZero(P);
    Copy(S,P.C(0,0));
  end Assign;

  function Scal(S: Scalar) return Taylor2 is
    P: Taylor2;
  begin
    SetZero(P);
    Copy(S,P.C(0,0));
    return P;
  end Scal;

  procedure SetZero(R: in Radius; P: in out Taylor2) is
  begin
    SetZero(P.C);
    P.F := FConst;
    P.R := R;
  end SetZero;

  procedure Coeff(I,J: in Natural; P: in Taylor2; S: in out Scalar) is
    F: constant Integer := TrueF(P);
    Deg: constant Integer := I+J;
    C: Poly2 renames P.C;
  begin
    if Deg <= F then
      if Deg <= PDeg then Copy(C(I,J),S); else SetZero(S); end if;
    else
      SetZero(S);
      if not STrunc then
        declare
          R: Radius renames P.R;
          E: Scalar;
        begin
          for D in F .. IMin(Deg-1,PDeg) loop      -- collect errors: degree >= D
            Div(R,S);
            for M in Imax(0,D-J) .. IMin(D,I) loop -- C(M,N) with M+N=D
              ModCenter(C(M,D-M),E);
              Add(E,S);
            end loop;
          end loop;
          if Deg <= PDeg then
            Div(R,S);
            Add(C(I,J),S);
          else
            Mult(Scal(R)**(PDeg-Deg),S);
          end if;
        end;
      end if;
    end if;
  end Coeff;

  function Coeff(I,J: Natural; P: Taylor2) return Scalar is
    S: Scalar;
  begin
    Coeff(I,J,P,S);
    return S;
  end Coeff;

  procedure AddCoeff(I,J: in Natural; S: in Scalar; P: in out Taylor2) is
    Deg: constant Integer := I+J;
  begin
    if not IsZero(S) then
      if Deg <= PDeg then
        Add(S,P.C(I,J));
        if P.F=FConst and then Deg>0 then P.F := PDeg1; end if;
      elsif not STrunc then
        declare
          C: Poly2 renames P.C;
          E: Scalar;
        begin
          Prod(Scal(P.R)**(PDeg-Deg),ToErr(S),E);
          if I<PDeg1Half then
            Add(E,P.C(I,PDeg-I));
          elsif J<PDeg1Half then
            Add(E,P.C(PDeg-J,J));
          else
            Mult(Half,E);
            Add(E,C(PDeg1Half,PDegHalf));
            Add(E,C(PDegHalf,PDeg1Half));
          end if;
          if P.F>PDeg then P.F := PDeg; end if;
        end;
      end if;
    end if;
  end AddCoeff;

  procedure Add(S: in Scalar; P: in out Taylor2) is
  begin
    Add(S,P.C(0,0));
  end Add;

  procedure Mult(S: in Scalar; P: in out Taylor2) is
  begin
    if IsZero(S) then SetZero(P); else Mult(S,P.C); end if;
  end Mult;

  procedure Prod(S: in Scalar; P1: in Taylor2; P2: in out Taylor2) is
  begin
    if IsZero(S) then
      SetZero(P2);
    else
      Prod(S,P1.C,P2.C);
      P2.F := P1.F;
      P2.R := P1.R;
    end if;
  end Prod;

  function "*"(S: Scalar; P: Taylor2) return Taylor2 is
    Q: Taylor2;
  begin
    Prod(S,P,Q);
    return Q;
  end "*";

  procedure AddProd(S: in Scalar; P1: in Taylor2; P2: in out Taylor2) is
  begin
    if not IsZero(S) then
      AddProd(S,P1.C,P2.C);
      if P2.F>P1.F then P2.F := P1.F; end if;
      if P2.R>P1.R then P2.R := P1.R; end if;
    end if;
  end AddProd;

  procedure SubProd(S: in Scalar; P1: in Taylor2; P2: in out Taylor2) is
  begin
    if not IsZero(S) then
      SubProd(S,P1.C,P2.C);
      if P2.F>P1.F then P2.F := P1.F; end if;
      if P2.R>P1.R then P2.R := P1.R; end if;
    end if;
  end SubProd;

  procedure Div(S: in Scalar; P: in out Taylor2) is
    Tmp: Scalar;
  begin
    Inv(S,Tmp);
    Mult(Tmp,P.C);
  end Div;

  procedure Val(P: in Taylor2; S1,S2: in Scalar; S3: in out Scalar) is
  begin
    if (P.F <= PDeg) and then (MaxNorm(S1)>P.R or else MaxNorm(S2)>P.R) then
      raise Undefined with Show0("Taylors2.Val error: argument too large");
    end if;
    Evaluate(P.C,S1,S2,S3);
  end Val;

  function Val(P: Taylor2; S1,S2: Scalar) return Scalar is
    S3: Scalar;
  begin
    Val(P,S1,S2,S3);
    return S3;
  end Val;

  procedure Val(P1: in Taylor2; S1,S2: in Scalar; P2: in out Taylor2) is
  begin
    SetZero(P2);
    Val(P1,S1,S2,P2.C(0,0));
  end Val;

  function NotSharp(C: Poly2; F,L: Integer) return Boolean is
  begin
    if not STrunc then
      for D in F .. L loop
        for I in 0 .. D loop
          if not IsSharp(C(I,D-I)) then return True; end if;
        end loop;
      end loop;
    end if;
    return False;
  end NotSharp;

  procedure Norm1(R: in Radius; P: in Taylor2; S: in out Scalar) is
    C: Poly2 renames P.C;
    Deg: constant Integer := QuasiDeg(C);
    Tmp: Scalar;
  begin
    if R>P.R and then NotSharp(C,P.F,Deg) then
      raise Domain_Error with Show0("Taylors2.Norm1: R>P.R");
    end if;
    SetZero(S);
    for D in reverse 0 .. Deg loop
      Mult(R,S,Tmp);
      for I in 0 .. D loop
        Norm(C(I,D-I),Tmp);
        Add(Tmp,S);
      end loop;
    end loop;
  end Norm1;

  procedure Norm1(P: in Taylor2; S: in out Scalar) is
  begin
    Norm1(P.R,P,S);
  end Norm1;

  function Norm1(R: Radius; P: Taylor2) return Scalar is
    S: Scalar;
  begin
    Norm1(R,P,S);
    return S;
  end Norm1;

  function Norm1(P: Taylor2) return Scalar is
    S: Scalar;
  begin
    Norm1(P.R,P,S);
    return S;
  end Norm1;

  procedure Monom(I,J: in Power; P: in out Taylor2) is
    D: constant Power := I+J;
  begin
    SetZero(P.C);
    Assign(1,P.C(I,J));
    if D=0 then P.F := FConst; else P.F := PDeg1; end if;
    P.R := Rad;
  end Monom;

  procedure HiUnitBall(I,J: in Power; P: in out Taylor2; R: in Radius := Rad) is
    D: constant Power := I+J;
  begin
    SetZero(P);
    P.C(I,J) :=  Scal(R)**(-D);
    ToErr(P.C(I,J));
    P.F := D;
    P.R := R;
  end HiUnitBall;

  procedure Compose(C: in SV.Polynom1; P1: in Taylor2; P2: in out Taylor2) is
    D: constant Integer := SV.EffLast(C);
    Tmp: LT_Pointer;
  begin
    SetZero(P2);
    Add(C(D),P2.C(0,0));
    if D>0 then
      Pool.Allocate(Tmp);
      for I in reverse 0 .. D-1 loop
        Mult(P1,P2,Tmp.all.Data);
        Add(C(I),P2.C(0,0));
      end loop;
      Pool.Recycle(Tmp);
    end if;
  end Compose;

  procedure Der1(R: in Radius; P: in out Taylor2) is
    C: Poly2 renames P.C;
  begin
    if P.F>PDeg then
      if P.F=FConst then
        SetZero(C);
      else
        Num_Der1(C);
        P.F := TrueF(P);
      end if;
      P.R := R;
      return;
    end if;
    if P.R <= R then
      raise Domain_Error with Show0("Taylors2.Der1: need P.R>R");
    end if;
    declare
      SM: constant Scalar := Inv(Log(Scal(P.R)/R));
      E1: constant Scalar := Exp(Scal(-1));
      RM: constant Flt := MaxNorm(SM);
      FM: constant Flt := MaxNorm(SM*E1/R);
      function DerFac(I: Positive) return Flt is
      begin
        if Flt(I)>RM then
          return ((R/P.R)**(I-1))/P.R;
        else
          return FM/Flt(I);
        end if;
      end DerFac;
      Fac: Flt;
      E: SV.Vector(Power);
    begin
      for J in P.F .. PDeg loop
        ModCenter(C(0,J),E(J));
      end loop;
      for I in Imax(P.F,1) .. PDeg loop
        for J in 0 .. PDeg-I loop
          ErrMult(DerFac(I),C(I,J));
        end loop;
      end loop;
      Num_Der1(C);
      Fac := DerFac(1)*P.R;
      for J in P.F .. PDeg loop
        AddProd(Fac,E(J),C(0,J));
      end loop;
    end;
    if P.F>0 then P.F := P.F-1; end if;
    P.R := R;
  end Der1;

  procedure Der2(R: in Radius; P: in out Taylor2) is
    C: Poly2 renames P.C;
  begin
    if P.F>PDeg then
      if P.F=FConst then
        SetZero(C);
      else
        Num_Der2(C);
        P.F := TrueF(P);
      end if;
      P.R := R;
      return;
    end if;
    if P.R <= R then
      raise Domain_Error with Show0("Taylors2.Der2: need P.R>R");
    end if;
    declare
      SM: constant Scalar := Inv(Log(Scal(P.R)/R));
      E1: constant Scalar := Exp(Scal(-1));
      RM: constant Flt := MaxNorm(SM);
      FM: constant Flt := MaxNorm(SM*E1/R);
      function DerFac(J: Positive) return Flt is
      begin
        if Flt(J)>RM then
          return ((R/P.R)**(J-1))/P.R;
        else
          return FM/Flt(J);
        end if;
      end DerFac;
      Fac: Flt;
      E: SV.Vector(Power);
    begin
      for I in P.F .. PDeg loop
        ModCenter(C(I,0),E(I));
      end loop;
      for J in Imax(P.F,1) .. PDeg loop
        for I in 0 .. PDeg-J loop
          ErrMult(DerFac(J),C(I,J));
        end loop;
      end loop;
      Num_Der2(C);
      Fac := DerFac(1)*P.R;
      for I in P.F .. PDeg loop
        AddProd(Fac,E(I),C(I,0));
      end loop;
    end;
    if P.F>0 then P.F := P.F-1; end if;
    P.R := R;
  end Der2;

  procedure ValDer1(P: in Taylor2; S1,S2: in Scalar; D1: in out Scalar) is
    Tmp: LT_Pointer;
  begin
    if P.F>PDeg then
      if P.F=FConst then
        SetZero(D1);
      else
        Num_ValDer1(P.C,S1,S2,D1);
      end if;
    else
      Pool.Allocate(Tmp);
      declare
        R: constant Radius := RMax(MaxNorm(S1),MaxNorm(S2));
        Q: Taylor2 renames Tmp.all.Data;
      begin
        Copy(P,Q);
        if R=Zero then
          LiftErrs(Q,1);
          Copy(Q.C(1,0),D1);
        else
          Der1(R,Q);
          Val(Q,S1,S2,D1);
        end if;
      end;
      Pool.Recycle(Tmp);
    end if;
  end ValDer1;

  procedure ValDer2(P: in Taylor2; S1,S2: in Scalar; D2: in out Scalar) is
    Tmp: LT_Pointer;
  begin
    if P.F>PDeg then
      if P.F=FConst then
        SetZero(D2);
      else
        Num_ValDer2(P.C,S1,S2,D2);
      end if;
    else
      Pool.Allocate(Tmp);
      declare
        R: constant Radius := RMax(MaxNorm(S1),MaxNorm(S2));
        Q: Taylor2 renames Tmp.all.Data;
      begin
        Copy(P,Q);
        if R=Zero then
          LiftErrs(Q,1);
          Copy(Q.C(0,1),D2);
        else
          Der2(R,Q);
          Val(Q,S1,S2,D2);
        end if;
      end;
      Pool.Recycle(Tmp);
    end if;
  end ValDer2;

  procedure ValDer(P: in Taylor2; S1,S2: in Scalar; D1,D2: in out Scalar) is
    Tmp: LT_Pointer;
  begin
    if P.F>PDeg then
      if P.F=FConst then
        SetZero(D1);
        SetZero(D2);
      else
        Num_ValDer1(P.C,S1,S2,D1);
        Num_ValDer2(P.C,S1,S2,D2);
      end if;
      return;
    end if;
    Pool.Allocate(Tmp);
    declare
      R: constant Radius := RMax(MaxNorm(S1),MaxNorm(S2));
      Q: Taylor2 renames Tmp.all.Data;
    begin
      Copy(P,Q);
      Der1(R,Q);
      Val(Q,S1,S2,D1);
      Copy(P,Q);
      Der2(R,Q);
      Val(Q,S1,S2,D2);
    end;
    Pool.Recycle(Tmp);
  end ValDer;

  procedure Num_MDer(P: in out Taylor2) is
  begin
    Num_MDer(P.C);
  end Num_MDer;

  procedure Inv_MDer(P: in out Taylor2; CheckConst: in Boolean := True) is
  begin
    Inv_MDer(P.C,CheckConst);
  end Inv_MDer;

  procedure NewQuot(S0,S1: in Scalar; Q: in out Flt; N: in out Integer) is
    use Flt_EF;
  begin
    if not (IsZero(S0) or else IsZero(S1)) then
      Q := Q+Log(MaxNorm(S0)/MaxNorm(S1));
      N := N+1;
    end if;
  end NewQuot;

  function GuessRho(P: Taylor2) return RadPair is
    use Flt_EF;
    DMin: constant Integer := PDeg/5;
    DMax: constant Integer := (4*PDeg)/5;
    C: Poly2 renames P.C;
    J: Integer;
    N1,N2: Integer := 0;
    Q1,Q2: Flt := Zero;
  begin
    for D in DMin .. DMax loop
      for I in 0 .. D loop
        J := D-I;
        NewQuot(C(I,J),C(I+1,J),Q1,N1);
        NewQuot(C(I,J),C(I,J+1),Q2,N2);
      end loop;
    end loop;
    return (Exp(Q1/Flt(N1)),Exp(Q2/Flt(N2)));
  end GuessRho;

  procedure HiCopy(P1: in Poly2; P2: in out Poly2) is
  begin
    for I in 0 .. PDeg loop
      Copy(P1(I,PDeg-I),P2(I,PDeg-I));
    end loop;
  end HiCopy;

  procedure HiErrMult(R: in Radius; P: in out Poly2) is
  begin
    for I in 0 .. PDeg loop
      ErrMult(R,P(I,PDeg-I));
    end loop;
  end HiErrMult;

  function HiContains(P1,P2: Poly2) return Logical is
    L: Logical := True;
  begin
    for I in 0 .. PDeg loop
      L := L and Contains(P1(I,PDeg-I),P2(I,PDeg-I));
      exit when L=False;
    end loop;
    return L;
  end HiContains;

  procedure Test_Some_Ops is
    R: constant Flt := 3.0/2.0;
    Fac: constant Flt := 1.0/4.0;
    X: constant Scalar := Scal(1.0/8.0);
    Y: constant Scalar := Scal(1.0/16.0);
    S0,S: Scalar;
    F0,F: Taylor2;
  begin
    SetZero(R,F0);
    Simple_Random(F0);
    Mult(Fac/MaxNorm(F0),F0);
    Add(1,F0);
    S0 := Val(F0,X,Y);              -- S0 := F0(X,Y);
    S := Inv(S0);
    F := Inv(F0);
    Show2("Inv(F(X,Y)) ",Val(F,X,Y),S);
    S := Sqrt(S0);
    F := Sqrt(F0);
    Show2("Sqrt(F(X,Y)) ",Val(F,X,Y),S);
    S := Root(5,S0);
    F := Root(5,F0);
    Show2("Root(5,F(X,Y)) ",Val(F,X,Y),S);
    S := Exp(S0);
    F := Exp(F0);
    Show2("Exp(F(X,Y)) ",Val(F,X,Y),S);
    S := Cos(S0);
    F := Cos(F0);
    Show2("Cos(F(X,Y)) ",Val(F,X,Y),S);
    S := Sin(S0);
    F := Sin(F0);
    Show2("Sin(F(X,Y)) ",Val(F,X,Y),S);
    S := Log(S0);
    F := Log(F0);
    Show2("Log(F(X,Y)) ",Val(F,X,Y),S);
    S := ArcCos(S0/Four);
    F := ArcCos(F0/Four);
    Show2("ArcCos((F(X,Y))/4) ",Val(F,X,Y),S);
  end Test_Some_Ops;

  function Num_MDer(P: in Taylor2) return Taylor2 is
    Q: Taylor2 := P;
  begin
    Num_MDer(Q.C);
    return Q;
  end Num_MDer;

  -------------------------------------------------
  -------------------------------------------------
  package TIP is new IPowers (Scalar => Taylor2);
  -------------------------------------------------
  -------------------------------------------------

  procedure IPower(I: in Integer; P1: in Taylor2; P2: in out Taylor2) renames TIP.IPower;

  function "**"(P: Taylor2; I: Integer) return Taylor2 is
    Q: Taylor2;
  begin
    IPower(I,P,Q);
    return Q;
  end "**";

  function Pi return Taylor2 is
  begin
    return Scal(Pi);
  end Pi;

  -----------------------------------------------
  --- stuff for Fun_Series ----------------------
  -----------------------------------------------

  procedure Inv_Split(P1: in Taylor2; P2: in out Taylor2; S: in out Scalar) is
  begin
    Inv(P1.C(0,0),S);
    Prod(S,P1,P2);
    SetZero(P2.C(0,0));
  end Inv_Split;

  procedure QPower_Split(Q: in Rational; P1: in Taylor2; P2: in out Taylor2; S: in out Scalar) is
  begin
    Copy(P1,P2);
    LiftErrs(P2,1);
    Inv(P2.C(0,0),S);
    Mult(S,P2);
    SetZero(P2.C(0,0));
    QPower(Q,P1.C(0,0),S);
  end QPower_Split;

  procedure Exp_Split(P1: in Taylor2; P2: in out Taylor2; S: in out Scalar) is
  begin
    Copy(P1,P2);
    LiftErrs(P2,1);
    Exp(P2.C(0,0),S);
    SetZero(P2.C(0,0));
  end Exp_Split;

  procedure CosSin_Split(P1: in Taylor2; P2: in out Taylor2; C,S: in out Scalar) is
  begin
    Copy(P1,P2);
    LiftErrs(P2,1);
    Cos(P2.C(0,0),C);
    Sin(P2.C(0,0),S);
    SetZero(P2.C(0,0));
  end CosSin_Split;

  procedure Log_Split(P1: in Taylor2; P2: in out Taylor2; S: in out Scalar) is
  begin
    Copy(P1,P2);
    LiftErrs(P2,1);
    Inv(P2.C(0,0),S);
    Mult(S,P2);
    SetZero(P2.C(0,0));
    Log(P1.C(0,0),S);
  end Log_Split;

  --------------------------------------------------------------------------------------
  --------------------------------------------------------------------------------------
  package TFS is new Fun_Series (Numeric => STrunc, Scalar => Scalar, Fun => Taylor2);
  --------------------------------------------------------------------------------------
  --------------------------------------------------------------------------------------

  procedure Inv(P1: in Taylor2; P2: in out Taylor2) is
  begin
    if P1.F=FConst then
      SetZero(P2);
      Inv(P1.C(0,0),P2.C(0,0));
    else
      declare
        Iter: constant Integer := Choose(STrunc,PDeg,PDeg+4);
      begin
        TFS.Series_Inv(P1,P2,Iter);
      end;
      if not STrunc then P2.F := IMin(P1.F,PDeg); end if;
    end if;
  end Inv;

  procedure Sqrt(P1: in Taylor2; P2: in out Taylor2) is
  begin
    if P1.F=FConst then
      SetZero(P2);
      Sqrt(P1.C(0,0),P2.C(0,0));
    else
      declare
        Iter: constant Integer := Choose(STrunc,PDeg,PDeg+4);
      begin
        TFS.Series_QPower(1/2,P1,P2,Iter);
      end;
      if not STrunc then P2.F := IMin(P1.F,PDeg); end if;
    end if;
  end Sqrt;

  procedure Root(K: in Positive; P1: in Taylor2; P2: in out Taylor2) is
  begin
    if K=1 then
      Copy(P1,P2);
    elsif P1.F=FConst then
      SetZero(P2);
      Root(K,P1.C(0,0),P2.C(0,0));
    else
      declare
        Iter: constant Integer := Choose(STrunc,PDeg,PDeg+4);
      begin
        TFS.Series_QPower(1/K,P1,P2,Iter);
      end;
      if not STrunc then P2.F := IMin(P1.F,PDeg); end if;
    end if;
  end Root;

  procedure QPower(Q: in Rational; P1: in Taylor2; P2: in out Taylor2) is
  begin
    if P1.F=FConst then
      SetZero(P2);
      QPower(Q,P1.C(0,0),P2.C(0,0));
    elsif IDen(Q)=1 then
      IPower(INum(Q),P1,P2);
    else
      declare
        Iter: constant Integer := Choose(STrunc,PDeg,PDeg+4);
      begin
        TFS.Series_QPower(Q,P1,P2,Iter);
      end;
      if not STrunc then P2.F := IMin(P1.F,PDeg); end if;
    end if;
  end QPower;

  procedure Exp(P1: in Taylor2; P2: in out Taylor2) is
  begin
    if P1.F=FConst then
      SetZero(P2);
      Exp(P1.C(0,0),P2.C(0,0));
    else
      declare
        Iter: constant Integer := Choose(STrunc,PDeg,PDeg+4);
      begin
        TFS.Series_Exp(P1,P2,Iter);
      end;
      if not STrunc then P2.F := IMin(P1.F,PDeg); end if;
    end if;
  end Exp;

  procedure CosSin(P: in Taylor2; PC,PS: in out Taylor2) is
  begin
    if P.F=FConst then
      SetZero(PC);
      Cos(P.C(0,0),PC.C(0,0));
      SetZero(PS);
      Sin(P.C(0,0),PS.C(0,0));
    else
      declare
        Iter: constant Integer := Choose(STrunc,PDeg/2,PDeg/2+4); -- uses P**2
      begin
        TFS.Series_CosSin(P,PC,PS,Iter);
      end;
      if not STrunc then PC.F := IMin(P.F,PDeg); PS.F := PC.F; end if;
    end if;
  end CosSin;

  procedure Log(P1: in Taylor2; P2: in out Taylor2) is
  begin
    if P1.F=FConst then
      SetZero(P2);
      Log(P1.C(0,0),P2.C(0,0));
    else
      declare
        Iter: constant Integer := Choose(STrunc,PDeg,PDeg+4);
      begin
        TFS.Series_Log(P1,P2,Iter);
      end;
      if not STrunc then P2.F := IMin(P1.F,PDeg); end if;
    end if;
  end Log;

  procedure ArcSin_Part(P1: in Taylor2; P2: in out Taylor2) is
    --- approx and modulo a constant
    Tmp: LT_Pointer;
  begin
    Pool.Allocate(Tmp);
    declare
      Q: Taylor2 renames Tmp.all.Data;
    begin
      Prod(P1,P1,Q);
      Center(Q);
      Neg(Q);
      Add(1,Q);
      QPower((-1)/2,Q,P2);   -- P2 := 1/Sqrt(1-P1*P1)
      Center(P2);
      Copy(P1,Q);
      Center(Q);
      Num_MDer(Q.C);
      Mult(Q,P2);            -- P2 := P1'/Sqrt(1-P1*P1)
    end;
    Pool.Recycle(Tmp);
    Center(P2);
    Inv_MDer(P2.C,False);    -- integrate P2
  end ArcSin_Part;

  --------------------------------------------------------------------------------------
  --------------------------------------------------------------------------------------
  package TN is new Newton(Numeric => STrunc, Scalar => Taylor2, Operator => Taylor2);
  --------------------------------------------------------------------------------------
  --------------------------------------------------------------------------------------

  procedure ArcCos(P1: in Taylor2; P2: in out Taylor2) is
  begin
    if P1.F=FConst then
      SetZero(P2);
      ArcCos(P1.C(0,0),P2.C(0,0));
      return;
    end if;
    ArcSin_Part(P1,P2);
    Neg(P2);
    ArcCos(P1.C(0,0),P2.C(0,0));
    Center(P2);
    if STrunc then return; end if;
    declare
      procedure F(X: in Taylor2; Y: in out Taylor2) is
      begin
        Cos(X,Y);
        Sub(P1,Y);
      end F;
      procedure DF(X: in Taylor2; Y: in out Taylor2) is
      begin
        Sin(X,Y);
        Neg(Y);
      end DF;
      procedure Improve is new TN.FindZero(F => F, DF => DF);
    begin
      Improve(P2,0);
    end;
    LiftErrs(P2,1);
    ArcCos(P1.C(0,0),P2.C(0,0));
  end ArcCos;

  procedure ArcSin(P1: in Taylor2; P2: in out Taylor2) is
  begin
    if P1.F=FConst then
      SetZero(P2);
      ArcSin(P1.C(0,0),P2.C(0,0));
      return;
    end if;
    ArcSin_Part(P1,P2);
    ArcSin(P1.C(0,0),P2.C(0,0));
    Center(P2);
    if STrunc then return; end if;
    declare
      procedure F(X: in Taylor2; Y: in out Taylor2) is
      begin
        Sin(X,Y);
        Sub(P1,Y);
      end F;
      procedure DF(X: in Taylor2; Y: in out Taylor2) is
      begin
        Cos(X,Y);
      end DF;
      procedure Improve is new TN.FindZero(F => F, DF => DF);
    begin
      Improve(P2,0);
    end;
    LiftErrs(P2,1);
    ArcSin(P1.C(0,0),P2.C(0,0));
  end ArcSin;

  ----------------------------------------------
  ----------------------------------------------
  package TR is new Roots (Scalar => Taylor2);
  ----------------------------------------------
  ----------------------------------------------

  procedure Roots2(B,C: in Taylor2; U1,U2,V: in out Taylor2) is
    Steps: constant Integer := 64;
  begin
    if Simple(B) and then Simple(C) then
      SetZero(U1);
      SetZero(U2);
      SetZero(V);
      Roots2(B.C(0,0),C.C(0,0),U1.C(0,0),U2.C(0,0),V.C(0,0));
    else
      TR.NewtonRoots2(B,C,U1,U2,V,Steps);
    end if;
  end Roots2;

  procedure Roots3(B,C,D: in Taylor2; U0,U1,U2,V: in out Taylor2) is
    Steps: constant Integer := IMax(PDeg,64);
  begin
    if Simple(B) and then Simple(C) and then Simple(D) then
      SetZero(U0);
      SetZero(U1);
      SetZero(U2);
      SetZero(V);
      Roots3(B.C(0,0),C.C(0,0),D.C(0,0),U0.C(0,0),U1.C(0,0),U2.C(0,0),V.C(0,0));
    else
      TR.NewtonRoots3(B,C,D,U0,U1,U2,V,Steps);
    end if;
  end Roots3;

  procedure Plot(Name: in String; G: in Taylor2; L: Positive := 50) is
    use Flt_IO;
    F: File_Type;
    X,Y,Z: Scalar;
  begin
    Create(F,Out_File,Name);
    for M in 0 .. L loop
      for N in 0 .. L loop
        Assign(G.R*Flt(2*M-L)/Flt(L),X);
        Assign(G.R*Flt(2*N-L)/Flt(L),Y);
        Val(G,X,Y,Z);
        Put(F,Approx(X),3,3,3);
        Put(F,Approx(Y),3,3,3);
        Put(F,Approx(Z),3,3,3);
        New_Line(F);
      end loop;
    end loop;
    Close(F);
  end Plot;

begin

  if FConst <= PDeg1 then
    raise Sorry with "Taylors2: increase FConst";
  end if;

end Taylors2;
