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

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

package body Taylors1 is

  use SM,SV,SP,SP_More;

  function NotConst(C: Poly1) return Boolean is
  begin
    for I in 1 .. PDeg loop
      if not IsZero(C(I)) then return True; end if;
    end loop;
    return False;
  end NotConst;

  procedure AddProd(S: in Scalar; C1: in Poly1; C2: in out Poly1; Tmp: in out Scalar) is
  begin
    if not IsZero(S) then
      for I in Power loop AddProd(S,C1(I),C2(I),Tmp); end loop;
    end if;
  end AddProd;

  procedure SubProd(S: in Scalar; C1: in Poly1; C2: in out Poly1; Tmp: in out Scalar) is
  begin
    if not IsZero(S) then
      for I in Power loop SubProd(S,C1(I),C2(I),Tmp); end loop;
    end if;
  end SubProd;

  --- basic

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

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

  function IsZero(P: Taylor1) return Boolean is
    C: Poly1 renames P.C;
  begin
    if P.F=FConst then
      if Check_Consistency and then NotConst(C) then raise Inconsistent_Data; end if;
      return IsZero(P.C(0));
    else
      return IsZero(P.C);
    end if;
  end IsZero;

  function "="(P1,P2: Taylor1) return Boolean is
    C1: Poly1 renames P1.C;
    C2: Poly1 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 Taylor1) is
  begin
    SetZero(P.C);
    P.F := FConst;
    P.R := Rad;
  end SetZero;

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

  procedure Swap(P1,P2: in out Taylor1) 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: Taylor1) return Boolean is
  begin
    return Center0(P.C);
  end Center0;

  function Contains0(P: Taylor1) return Logical is
    --- very crude
    L: Logical := Contains0(P.C(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: Taylor1) return Logical is
    --- very crude
    C1: Poly1 renames P1.C;
    C2: Poly1 renames P2.C;
    L: Logical := Contains(C1(0),C2(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 Taylor1) is
  begin
    SetZero(P.C);
    if STrunc then
      P.F := FConst;
    else
      BallAt0(R,P.C(0));
      P.F := 0;
    end if;
    P.R := Rad;
  end BallAt0;

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

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

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

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

  procedure Center(P: in out Taylor1) is
    C: Poly1 renames P.C;
  begin
    if Not_STrunc then
      if P.F=FConst then
        if Check_Consistency and then NotConst(C) then raise Inconsistent_Data; end if;
        Center(C(0));
      else
        Center(C);
        if EffDeg0(C) then
          P.F := FConst;
        else
          P.F := PDeg1;
        end if;
      end if;
    end if;
  end Center;

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

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

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

  procedure ModCenter(P1: in Taylor1; P2: in out Taylor1) 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: Taylor1) return Taylor1 is
    Q: Taylor1;
  begin
    ModCenter(P,Q);
    return Q;
  end ModCenter;

  procedure ErrMult(R: in Radius; P: in out Taylor1) is
    C: Poly1 renames P.C;
  begin
    if Not_STrunc then for I in Power loop ErrMult(R,C(I)); end loop; end if;
  end ErrMult;

  procedure Union(P1: in Taylor1; P2: in out Taylor1) 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: Taylor1) return Taylor1 is
    Q: Taylor1;
  begin
    Copy(P2,Q);
    Union(P1,Q);
    return Q;
  end Union;

  procedure Intersection(P1: in Taylor1; P2: in out Taylor1; Empty: out Logical) is
    --- far from optimal
    C1: Poly1 renames P1.C;
    C2: Poly1 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
      TPool.Allocate(Tmp);
      Copy(P2,Tmp.all.Data);
      Copy(P1,P2);
      Intersection(Tmp.all.Data,P2,Empty);
      TPool.Recycle(Tmp);
      return;
    end if;
    if P2.R<P1.R then P2.R := P1.R; end if;
    Empty := False;
    for I in 0 .. IMin(F1,PDeg) loop
      Intersection(C1(I),C2(I),E);
      Empty := Empty or E;
    end loop;
    if F1<PDeg then Empty := Empty or Uncertain; end if;
  end Intersection;

  function Interval(C: Flt; R: Radius) return Taylor1 is
     P: Taylor1;
  begin
     SetZero(P);
     P.R := Rad;
     if R=Zero then
        P.F := FConst;
        Assign(C,P.C(0));
     elsif PDeg=0 then
        P.F := FConst;
        Enclose(-(R-C),C+R,P.C(0));
     else
        P.F := PDeg1;
        Assign(C,P.C(0));
        Assign(R,P.C(1));
     end if;
     return P;
  end Interval;

  --- order

  function Sign(P: Taylor1) return Integer is
    Prec: constant Radius := 1.0E-2;
  begin
    if not IsZero(P) then
      if not IsReal(P) then
        raise Undefined with Show0("Taylor1.Sign error: argument non-real");
      end if;
      if Inf(P,Prec)>Zero then return  1; end if;
      if Sup(P,Prec)<Zero then return -1; end if;
      raise Not_Certain;
    end if;
    return 0;
  end Sign;

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

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

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

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

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

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

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

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

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

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

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

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

  function QuasiDeg(P: Taylor1) return Natural is
  begin
    if P.F=FConst then
      if Check_Consistency and then NotConst(P.C) then raise Inconsistent_Data; end if;
      return 0;
    else
      return EffDeg(P.C);
    end if;
  end QuasiDeg;

  function GuessMaxAbs(P: Taylor1; D: Natural) return Radius is
    KMax: constant Integer := 1+D/2;
    C: Poly1 renames P.C;
    R: Flt renames P.R;
    X,N: Flt;
    Y: Scalar;
  begin
    N := MaxNorm(C(0));
    for K in -KMax .. KMax loop
      if K /= 0 then
        X := R*Flt(K)/Flt(KMax);
        Copy(C(D),Y);                   -- inlined Evaluate
        for I in reverse 0 .. D-1 loop  -- inlined Evaluate
          Mult(X,Y);                    -- inlined Evaluate
          Add(C(I),Y);                  -- inlined Evaluate
        end loop;                       -- inlined Evaluate Y := p(X)
        N := RMax(MaxNorm(Y),N);
      end if;
    end loop;
    return N;
  end GuessMaxAbs;

  function Sharp_Sup(P: Taylor1; Prec: Radius) return Flt is
    D: constant Natural := QuasiDeg(P);
    C: Poly1 renames P.C;
    R: Radius renames P.R;
    V: Flt := Sup(C(0));
  begin
    if D<2 then
      if D=0 then return V; end if;
      return RMax(V+R*Sup(C(1)),V-R*Inf(C(1)));
    end if;
    declare
      Eps: constant Radius := Prec*GuessMaxAbs(P,D);
      W: Flt;
      Y,Tmp: Scalar;

      procedure Val1(X: in Flt) is     -- Y := p(X)
      begin
        Copy(C(D),Y);
        for I in reverse 0 .. D-1 loop
          Mult(X,Y);
          Add(C(I),Y);
        end loop;
      end Val1;

      procedure Val2(X1,X2: in Flt) is  -- Y := p([X1,X2])
      begin
        Enclose(X1,X2,Tmp);
        Copy(C(D),Y);
        for I in reverse 0 .. D-1 loop
          Mult(Tmp,Y);
          Add(C(I),Y);
        end loop;
      end Val2;

      procedure FindSup(L: in Integer; X1,X2: in Flt) is
      begin
        if STrunc then
          Val1(X2);    -- Y := p(X2)
          W := Sup(Y);
          if W <= V then return; end if;
          Swap(Y,Tmp);
          Val1(X1);    -- Y := p(X1)
          Sub(Y,Tmp);
        else
          Val2(X1,X2); -- Y := p([X1,X2])
          W := Sup(Y);
          if W <= V then return; end if;
          ModCenter(Y,Tmp);
        end if;
        if MaxNorm(Tmp)<Eps then
          V := W;
        elsif L=0 then
          raise Sorry with "Taylors1.Sup: getting tired";
        else
          FindSup(L-1,X1,Half*(X1+X2)); -- half-interval by X1
          FindSup(L-1,Half*(X1+X2),X2); -- half-interval by X2
        end if;
      end FindSup;

      L: constant Integer := 128; --- maximal recursion level
      X0,X: Flt := Zero;
    begin
      for K in 1 .. D loop
        X0 := X;
        if K=D then X := R; else X := R*Flt(K)/Flt(D); end if;
        FindSup(L, X0, X);
        FindSup(L,-X0,-X);
      end loop;
      return V;
    end;
  end Sharp_Sup;

  function Sup(P: Taylor1; Prec: Radius) return Flt is
    --- assumes Flt operations are rounded up
  begin
    if IsSharp(P) then
      return Sharp_Sup(P,Prec);
    else
      declare
        Tmp: LT_Pointer;
        V: Flt;
      begin
        TPool.Allocate(Tmp);
        Center(P,Tmp.all.Data);
        V := Sup(Tmp.all.Data,Prec);
        ModCenter(P,Tmp.all.Data);
        V := V+MaxNorm(Tmp.all.Data);
        TPool.Recycle(Tmp);
        return V;
      end;
    end if;
  end Sup;

  function Sup(P: Taylor1) return Flt is
  begin
    return Sup(P,Psi);
  end Sup;

  function Inf(P: Taylor1; Prec: Radius) return Flt is
    R: Flt;
    Tmp: LT_Pointer;
  begin
    TPool.Allocate(Tmp);
    Neg(P,Tmp.all.Data);
    R := Sup(Tmp.all.Data,Prec);
    TPool.Recycle(Tmp);
    return -R;
  end Inf;

  function Inf(P: Taylor1) return Flt is
  begin
    return Inf(P,Psi);
  end Inf;

  --- addition and multiplication etc.

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

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

  function "-"(P: Taylor1) return Taylor1 is
    Q: Taylor1;
  begin
    Neg(P.C,Q.C);
    Q.F := P.F;
    Q.R := P.R;
    return Q;
  end "-";

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

  procedure Add(P1: in Taylor1; P2: in out Taylor1) is
  begin
    if P2.R>P1.R then
      if P1.F=FConst then
        if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
        Add(P1.C(0),P2.C(0));
        return;
      end if;
      if P2.F <= PDeg then P2.R := P1.R; end if;
    end if;
    if P2.F>P1.F then P2.F := P1.F; end if;
    Add(P1.C,P2.C);
  end Add;

  procedure Sum(P1,P2: in Taylor1; P3: in out Taylor1) is
  begin
    if P1.F=FConst then
      if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
      Copy(P2,P3);
      Add(P1.C(0),P3.C(0));
    elsif P2.F=FConst then
      if Check_Consistency and then NotConst(P2.C) then raise Inconsistent_Data; end if;
      Copy(P1,P3);
      Add(P2.C(0),P3.C(0));
    else
      Sum(P1.C,P2.C,P3.C);
      P3.F := IMin(P1.F,P2.F);
      P3.R := RMin(P1.R,P2.R);
    end if;
  end Sum;

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

  procedure Sub(P1: in Taylor1; P2: in out Taylor1) is
  begin
    if P2.R>P1.R then
      if P1.F=FConst then
        if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
        Sub(P1.C(0),P2.C(0));
        return;
      end if;
      if P2.F <= PDeg then P2.R := P1.R; end if;
    end if;
    if P2.F>P1.F then P2.F := P1.F; end if;
    Sub(P1.C,P2.C);
  end Sub;

  procedure Diff(P1,P2: in Taylor1; P3: in out Taylor1) is
  begin
    if P1.F=FConst then
      if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
      Neg(P2,P3);
      Add(P1.C(0),P3.C(0));
    elsif P2.F=FConst then
      if Check_Consistency and then NotConst(P2.C) then raise Inconsistent_Data; end if;
      Copy(P1,P3);
      Sub(P2.C(0),P3.C(0));
    else
      Diff(P1.C,P2.C,P3.C);
      P3.F := IMin(P1.F,P2.F);
      P3.R := RMin(P1.R,P2.R);
    end if;
  end Diff;

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

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

  procedure Mult(R: in Flt; P: in out Taylor1; Tmp: in out Scalar) is
  begin
    if R=Zero then
      SetZero(P);
    else
      Assign(R,Tmp);
      Mult(Tmp,P);
    end if;
  end Mult;

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

  procedure Prod(R: in Flt; P1: in Taylor1; P2: in out Taylor1) is
  begin
    if R=Zero then
      SetZero(P2);
    else
      declare
        STmp: Scalar;
      begin
        Assign(R,STmp);
        Prod(STmp,P1,P2);
      end;
    end if;
  end Prod;

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

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

  procedure AddProd(R: in Flt; P1: in Taylor1; P2: in out Taylor1) is
  begin
    if R /= Zero then
      declare
        STmp: Scalar;
      begin
        Assign(R,STmp);
        AddProd(STmp,P1,P2);
      end;
    end if;
  end AddProd;

  procedure Mult(Q: in Rational; P: in out Taylor1) is
  begin
    if LNum(Q)=0 then
      SetZero(P);
    else
      declare
        Tmp: Scalar;
      begin
        Mult(Flt(LNum(Q)),P,Tmp);
        Div(Flt(LDen(Q)),P,Tmp);
      end;
    end if;
  end Mult;

  function Prod_Err(C1,C2: Poly1; D1,D2: Natural; R: Radius) return Radius is
    --- uses FBalls
    D3: constant Natural := D1+D2;
    E: Radius := Zero;
    S: FBall;
    L1,L2: LB_Pointer;
  begin
    BPool.Allocate(L1);
    BPool.Allocate(L2);
    declare
      P1: FPoly1 renames L1.all.Data;
      P2: FPoly1 renames L2.all.Data;
    begin
      for I in 0 .. D1 loop
        Enclose(Inf(C1(I)),Sup(C1(I)),P1(I));
      end loop;
      for I in 0 .. D2 loop
        Enclose(Inf(C2(I)),Sup(C2(I)),P2(I));
      end loop;
      for D in reverse PDeg1 .. D3 loop
        S := (Zero,Zero);
        for I in D-D2 .. D1 loop
          AddProd(P1(I),P2(D-I),S);
        end loop;
        E := R*(E+MaxNorm(S));
      end loop;
    end;
    BPool.Recycle(L1);
    BPool.Recycle(L2);
    return E;
  end Prod_Err;

  procedure Add_Prod(P,Q: in Taylor1; H: in out Taylor1) is
    PC: Poly1 renames P.C;
    QC: Poly1 renames Q.C;
    HC: Poly1 renames H.C;
    P0: Scalar renames PC(0);
    Q0: Scalar renames QC(0);
  begin
    if P.F=FConst or else PDeg=0 then AddProd(P0,Q,H); return; end if;
    if Q.F=FConst or else PDeg=0 then AddProd(Q0,P,H); return; end if;
    if H.F>P.F then H.F := P.F; end if;
    if H.F>Q.F then H.F := Q.F; end if;
    if H.R>P.R then H.R := P.R; end if;
    if H.R>Q.R then H.R := Q.R; end if;
    if PDeg<3 then
      declare
        P1: Scalar renames PC(1);
        Q1: Scalar renames QC(1);
        H0: Scalar renames HC(0);
        H1: Scalar renames HC(1);
        Tmp: Scalar;
      begin
        AddProd(P0,Q0,H0,Tmp);
        AddProd(P0,Q1,H1,Tmp);
        AddProd(P1,Q0,H1,Tmp);
        if PDeg=1 then                      -- case PDeg=1
          if Not_STrunc then
            declare
              E: Scalar;
            begin
              Prod(P1,Q1,E);
              ToErr(E);
              ErrMult(Q.R,E);
              Add(E,H1);
            end;
          end if;
          return;
        end if;
        pragma Warnings (Off);
        declare                              -- case PDeg=2
          P2: Scalar renames PC(2);
          Q2: Scalar renames QC(2);
          H2: Scalar renames HC(2);
        begin
          AddProd(P0,Q2,H2,Tmp);
          AddProd(P1,Q1,H2,Tmp);
          AddProd(P2,Q0,H2,Tmp);
          if Not_STrunc then
            declare
              E: Scalar;
            begin
              Prod(P2,Q2,E);
              ToErr(E);
              ErrMult(Q.R,E);
              AddProd(P1,Q2,E,Tmp);
              AddProd(P2,Q1,E,Tmp);
              ToErr(E);
              ErrMult(Q.R,E);
              Add(E,H2);
            end;
          end if;
        end;
      end;
      return;
    end if;
    pragma Warnings (On);
    declare
      DP: constant Integer := EffDeg(PC);  -- case PDeg>2
      DQ: constant Integer := EffDeg(QC);
      DH: constant Integer := DP+DQ;
    begin
      if DP=0 then AddProd(P0,Q,H); return; end if;
      if DQ=0 then AddProd(Q0,P,H); return; end if;
      declare
        Tmp: Scalar;
      begin
        for D in 0 .. IMin(DH,PDeg) loop
          declare
            HCD: Scalar renames HC(D);
          begin
            for I in 0 .. D loop
              AddProd(PC(I),QC(D-I),HCD,Tmp);
            end loop;
          end;
        end loop;
        if STrunc or else DH <= PDeg then return; end if;
        if Use_FBalls then
          BallAt0(Prod_Err(PC,QC,DP,DQ,H.R),Tmp); -- higher order via FBalls
          Add(Tmp,HC(PDeg));
          return;
        end if;
        declare
          S,E: Scalar;
        begin                                    -- higher order
          SetZero(E);
          for D in reverse PDeg1 .. DH loop
            SetZero(S);
            for I in D-DQ .. DP loop
              AddProd(PC(I),QC(D-I),S,Tmp);
            end loop;
            ToErr(S);
            Add(S,E);
            Mult(H.R,E,Tmp);
          end loop;
          Add(E,HC(PDeg));
        end;
      end;
    end;
  end Add_Prod;

  procedure Mult(P1: in Taylor1; P2,Tmp: in out Taylor1) is
  begin
    if P1.F=FConst then
      if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
      Mult(P1.C(0),P2);
    elsif P2.F=FConst then
      if Check_Consistency and then NotConst(P2.C) then raise Inconsistent_Data; end if;
      declare
        STmp: Scalar renames Tmp.C(0);
      begin
        Copy(P2.C(0),STmp);
        Prod(STmp,P1,P2);
      end;
    else
      Copy(P2,Tmp);
      SetZero(P2.C);
      Add_Prod(P1,Tmp,P2);
    end if;
  end Mult;

  procedure Mult(P1: in Taylor1; P2: in out Taylor1) is
    PTmp: LT_Pointer;
  begin
    if P1.F=FConst then
      if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
      Mult(P1.C(0),P2);
    elsif P2.F=FConst then
      if Check_Consistency and then NotConst(P2.C) then raise Inconsistent_Data; end if;
      declare
        STmp: Scalar;
      begin
        Copy(P2.C(0),STmp);
        Prod(STmp,P1,P2);
      end;
    else
      TPool.Allocate(PTmp);
      Copy(P2,PTmp.all.Data);
      SetZero(P2.C);
      Add_Prod(P1,PTmp.all.Data,P2);
      TPool.Recycle(PTmp);
    end if;
  end Mult;

  procedure Prod(P1,P2: in Taylor1; P3: in out Taylor1) is
  begin
    if P1.F=FConst then
      if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
      Prod(P1.C(0),P2,P3);
    elsif P2.F=FConst then
      if Check_Consistency and then NotConst(P2.C) then raise Inconsistent_Data; end if;
      Prod(P2.C(0),P1,P3);
    else
      SetZero(P1.R,P3);
      Add_Prod(P1,P2,P3);
    end if;
  end Prod;

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

  procedure AddProd(P1,P2: in Taylor1; P3: in out Taylor1) is
  begin
    if P1.F=FConst then
      if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
      AddProd(P1.C(0),P2,P3);
    elsif P2.F=FConst then
      if Check_Consistency and then NotConst(P2.C) then raise Inconsistent_Data; end if;
      AddProd(P2.C(0),P1,P3);
    else
      Add_Prod(P1,P2,P3);
    end if;
  end AddProd;

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

  procedure SubProd(P1,P2: in Taylor1; P3: in out Taylor1) is
  begin
    if P1.F=FConst then
      if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
      SubProd(P1.C(0),P2,P3);
    elsif P2.F=FConst then
      if Check_Consistency and then NotConst(P2.C) then raise Inconsistent_Data; end if;
      SubProd(P2.C(0),P1,P3);
    else
      Neg(P3);
      Add_Prod(P1,P2,P3);
      Neg(P3);
    end if;
  end SubProd;

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

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

  procedure Div(R: in Flt; P: in out Taylor1; Tmp: in out Scalar) is
  begin
    Quot(SOne,R,Tmp);
    Mult(Tmp,P);
  end Div;

  procedure Div(R: in Flt; P: in out Taylor1) is
    STmp: Scalar;
  begin
    Quot(SOne,R,STmp);
    Mult(STmp,P);
  end Div;

  procedure Quot(P1: in Taylor1; R: in Flt; P2: in out Taylor1) is
    STmp: Scalar;
  begin
    Quot(SOne,R,STmp);
    Prod(STmp,P1,P2);
  end Quot;

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

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

  procedure Quot(P1,P2: in Taylor1; P3: in out Taylor1) is
  begin
    Inv(P2,P3);
    Mult(P1,P3);
  end Quot;

  function "/"(P1,P2: Taylor1) return Taylor1 is
    P: Taylor1;
  begin
    Inv(P2,P);
    Mult(P1,P);
    return P;
  end "/";

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

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

  --- functions

  function IsReal(P: Taylor1) return Boolean is
  begin
    if TInfo.IsComplex then
      if P.F=FConst then
        if Check_Consistency and then NotConst(P.C) then raise Inconsistent_Data; end if;
        return IsReal(P.C(0));
      else
        return IsReal(P.C);
      end if;
    else
      return True;
    end if;
  end IsReal;

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

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

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

  procedure Real_Part(P: in out Taylor1) 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 Taylor1) 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 Taylor1) is
    C: Poly1 renames P.C;
  begin
    Eval0(C(0));
    for I in 1 .. PDeg loop SetZero(C(I)); end loop;
    P.F := FConst;
  end Eval0;

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

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

  function MaxNorm(P: Taylor1) return Radius is
    --- assumes Flt operations are rounded up
    C: Poly1 renames P.C;
    R: Radius renames P.R;
    N: Flt := Zero;
  begin
    for I in reverse 0 .. EffDeg(C) loop
      N := R*N+MaxNorm(C(I));
    end loop;
    return N;
  end MaxNorm;

  procedure Sqr(P: in out Taylor1) is
    Tmp: LT_Pointer;
  begin
    if P.F=FConst then
      if Check_Consistency and then NotConst(P.C) then raise Inconsistent_Data; end if;
      Sqr(P.C(0));
    else
      TPool.Allocate(Tmp);
      Copy(P,Tmp.all.Data);
      Mult(Tmp.all.Data,P);
      TPool.Recycle(Tmp);
    end if;
  end Sqr;

  function Sqr(P: Taylor1) return Taylor1 is
    Q: Taylor1;
  begin
    Copy(P,Q);
    Mult(P,Q);
    return Q;
  end Sqr;

  procedure Sqr(P1: in Taylor1; P2: in out Taylor1) is
  begin
    SetZero(P1.R,P2);
    if P1.F=FConst then
      if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
      Sqr(P1.C(0),P2.C(0));
    else
      Add_Prod(P1,P1,P2);
    end if;
  end Sqr;

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

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

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

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

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

  procedure Cos(P1: in Taylor1; P2: in out Taylor1) is
    Tmp: LT_Pointer;
  begin
    if P1.F=FConst then
      if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
      SetZero(P2);
      Cos(P1.C(0),P2.C(0));
    else
      TPool.Allocate(Tmp);
      CosSin(P1,P2,Tmp.all.Data);
      TPool.Recycle(Tmp);
    end if;
  end Cos;

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

  procedure Sin(P1: in Taylor1; P2: in out Taylor1) is
    Tmp: LT_Pointer;
  begin
    if P1.F=FConst then
      if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
      SetZero(P2);
      Sin(P1.C(0),P2.C(0));
    else
      TPool.Allocate(Tmp);
      CosSin(P1,Tmp.all.Data,P2);
      TPool.Recycle(Tmp);
    end if;
  end Sin;

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

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

  procedure Simple_Random(P: in out Taylor1) is
    R: constant Radius := P.R*Flt(9)/Flt(8); -- ad hoc
    C: Poly1 renames P.C;
    PowR: Flt := One;
  begin
    Simple_Random(C(0));
    for I in 1 .. PDeg loop
      PowR := PowR/R;
      Simple_Random(C(I));
      Mult(PowR,C(I));
    end loop;
    if PDeg>0 then P.F := PDeg1; end if;
  end Simple_Random;

  --- conversion and i/o

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

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

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

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

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

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

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

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

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

  procedure Show1(N: in String; P: in Taylor1; 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 Taylor1; 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 Taylor1; Decimal: in Boolean := False) is
    C: Poly1 renames P.C;
    Deg: constant Integer := EffDeg(C);
    PF:  constant Integer := IMin(TrueF(P),Deg+1);
  begin
    Put(F,PF);
    Put(F,P.R,Decimal);
    Put(F,Deg);
    for I in 0 .. Deg loop
      Put(F,C(I),Decimal);
    end loop;
  end Put;

  procedure Get(F: in File_Type; P: in out Taylor1; 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 I in 0 .. Deg loop
      Get(F,S,Decimal);
      AddCoeff(I,S,P);
    end loop;
    P.F := TrueF(P);
  end Get;

  procedure Write(FileName: in String; P: in Taylor1; 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 Taylor1; Decimal: in Boolean := False) is
    F: File_Type;
  begin
    if Verbosity>1 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 Taylor1 is
    P: Taylor1;
  begin
    Read(FileName,P,Decimal);
    return P;
  end Read;

  --- misc

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

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

  procedure Free_Cache(Dummy: in Taylor1) is
    pragma Unreferenced(Dummy);
  begin
    null;
  end Free_Cache;

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

  procedure Show_Param_T1(N: in String := "") is
  begin
    if STrunc then
      Show0(N & "T1 Scalar type is Numeric");
    else
      Show0(N & "T1 Scalar type is Ball");
    end if;
    Show1(N & "T1 Precision ",Get_Precision(SOne/Three));
    Show1(N & "T1 PDeg = ",PDeg);
  end Show_Param_T1;

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

  procedure SetRho(R: in Radius; P: in out Taylor1; 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("Taylors1.SetRho: cannot increase domain");
    end if;
    P.R := R;
  end SetRho;

  function TrueF(P: Taylor1) return Natural is
    --- the value FConst means that P is constant
    --- the value PDeg1 means that P is a polynomial
    C: Poly1 renames P.C;
  begin
    if P.F=FConst then
      if Check_Consistency and then NotConst(C) then raise Inconsistent_Data; end if;
      return FConst;
    end if;
    declare
      Deg: constant Integer := EffDeg(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)) then return FConst; end if;
        return 0;
      end if;
      for I in P.F .. Deg loop
        if not IsSharp(C(I)) then return I; end if;
      end loop;
      return PDeg1;
    end;
  end TrueF;

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

  procedure AdjustF(P: in out Taylor1; D: out Natural) is
    F: Natural renames P.F;
    C: Poly1 renames P.C;
  begin
    D := EffDeg(C);
    if D=0 then
      if STrunc or else F>0 or else IsSharp(C(0)) then
        F := FConst;
      else
        F := 0;
      end if;
    elsif F>D then
      if F=FConst then raise Inconsistent_Data; end if;
      F := PDeg1;
    elsif STrunc then
      F := PDeg1;
    else
      while IsSharp(C(F)) loop
        F := F+1;
        exit when F>D;
      end loop;
    end if;
  end AdjustF;

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

  procedure LiftErrs(P: in out Taylor1; F: in Power := 0) is
  begin
    P.F := TrueF(P);
    if F>P.F then
      declare
        R: Radius renames P.R;
        C: Poly1 renames P.C;
        E: Scalar;
      begin
        for K in P.F .. F-1 loop
          ModCenter(C(K),E);
          Div(R,E);
          Add(E,C(K+1));
        end loop;
      end;
      P.F := F;
    end if;
  end LiftErrs;

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

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

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

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

  procedure SetOne(R: in Radius; P: in out Taylor1) is
    C: Poly1 renames P.C;
  begin
    SetZero(C);
    Copy(SOne,C(0));
    P.F := FConst;
    P.R := R;
  end SetOne;

  procedure Affine(C0,C1: in Scalar; P: in out Taylor1) is
    C: Poly1 renames P.C;
  begin
    SetZero(C);
    Copy(C0,C(0));
    Copy(C1,C(1));
    P.F := PDeg1;
    P.R := Rad;
  end Affine;

  procedure Component(I: in Natural; P: in Taylor1; S: in out Scalar) is
    C: Poly1 renames P.C;
  begin
    if I=0 then
      Copy(C(0),S);
    elsif P.F=FConst then
      if Check_Consistency and then NotConst(P.C) then raise Inconsistent_Data; end if;
      SetZero(S);
    elsif I>PDEg then
      SetZero(S);
    else
      Copy(C(I),S);
    end if;
  end Component;

  procedure Coeff(I: in Natural; P: in Taylor1; S: in out Scalar) is
    F: constant Integer := TrueF(P);
    C: Poly1 renames P.C;
  begin
    if I <= F then
      if I <= PDeg then Copy(C(I),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(I-1,PDeg) loop
            Div(R,S);
            ModCenter(C(D),E);
            Add(E,S);
          end loop;
          if I <= PDeg then
            Div(R,S);
            Add(C(I),S);
          else
            Mult(Scal(R)**(PDeg-I),S);
          end if;
        end;
      end if;
    end if;
  end Coeff;

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

  procedure AddCoeff(I: in Natural; S: in Scalar; P: in out Taylor1) is
    --- assuming Flt ops are rounded up
  begin
    if not IsZero(S) then
      if I <= PDeg then
        Add_Coeff(I,S,P);
      elsif Not_STrunc then
        AddProd(P.R**(I-PDeg),ToErr(S),P.C(PDeg));
        if P.F>PDeg then P.F := PDeg; end if;
      end if;
    end if;
  end AddCoeff;

  procedure Add_Coeff(I: in Power; S: in Scalar; P: in out Taylor1) is
    --- assuming Flt ops are rounded up
  begin
    Add(S,P.C(I));
    if P.F=FConst then
      if Check_Consistency and then NotConst(P.C) then raise Inconsistent_Data; end if;
      if I>0 then P.F := PDeg1; end if;
    end if;
  end Add_Coeff;

  procedure Sub_Coeff(I: in Power; S: in Scalar; P: in out Taylor1) is
    --- assuming Flt ops are rounded up
  begin
    Sub(S,P.C(I));
    if P.F=FConst then
      if Check_Consistency and then NotConst(P.C) then raise Inconsistent_Data; end if;
      if I>0 then P.F := PDeg1; end if;
    end if;
  end Sub_Coeff;

  procedure Extract_Coeff(I: in Power; P: in out Taylor1; S: in out Scalar) is
    C: Poly1 renames P.C;
  begin
    if STrunc then
      SetZero(S);
      Swap(C(I),S);
    elsif I=PDeg then
      LiftErrs(P,I);
      Copy(C(I),S);
      ModCenter(C(I));
    else
      LiftErrs(P,I+1);
      SetZero(S);
      Swap(C(I),S);
    end if;
  end Extract_Coeff;

  procedure Add_Ball(I: in Power; R: in Flt; P: in out Taylor1) is
  begin
    if not (STrunc or else R=Zero) then
      declare
        B: Scalar;
      begin
        BallAt0(R,B);
        Add(B,P.C(I));
        if P.F>I then P.F := I; end if;
      end;
    end if;
  end Add_Ball;

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

  procedure Mult(S: in Scalar; P: in out Taylor1) is
  begin
    if IsZero(S) then
      SetZero(P.C);
      P.F := FConst;
    else
      declare
        C: Poly1 renames P.C;
      begin
        for I in Power loop Mult(S,C(I)); end loop;
      end;
    end if;
  end Mult;

  procedure Prod(S: in Scalar; P1: in Taylor1; P2: in out Taylor1) is
  begin
    if IsZero(S) then
      SetZero(P2.C);
      P2.F := FConst;
    else
      declare
        C1: Poly1 renames P1.C;
        C2: Poly1 renames P2.C;
      begin
        for I in Power loop Prod(S,C1(I),C2(I)); end loop;
      end;
      P2.F := P1.F;
    end if;
    P2.R := P1.R;
  end Prod;

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

  procedure AddProd(S: in Scalar; P1: in Taylor1; P2: in out Taylor1) is
  begin
    if not IsZero(S) then
      declare
        Tmp: Scalar;
      begin
        AddProd(S,P1.C,P2.C,Tmp);
      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 SubProd(S: in Scalar; P1: in Taylor1; P2: in out Taylor1) is
  begin
    if not IsZero(S) then
      declare
        Tmp: Scalar;
      begin
        SubProd(S,P1.C,P2.C,Tmp);
      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 SubProd;

  procedure Div(S: in Scalar; P: in out Taylor1) is
    C: Poly1 renames P.C;
    Si: Scalar;
  begin
    Inv(S,Si);
    for I in Power loop Mult(Si,C(I)); end loop;
  end Div;

  procedure Val(P: in Taylor1; S1: in Scalar; S2: in out Scalar) is
    D: Natural;
    C: Poly1 renames P.C;
  begin
    if IsZero(S1) then
      Copy(C(0),S2);
    else
      if (P.F <= PDeg) and then MaxNorm(S1)>P.R then
        Show2("",MaxNorm(S1),P.R);
        raise Undefined with "Taylor1s.Val error: argument too large";
      end if;
      D := EffLast(C);
      Copy(C(D),S2);
      for I in reverse 0 .. D-1 loop
        Mult(S1,S2);
        Add(C(I),S2);
      end loop;
    end if;
  end Val;

  function Val(P: Taylor1; S: Scalar) return Scalar is
    V: Scalar;
  begin
    Val(P,S,V);
    return V;
  end Val;

  procedure Val(P1: in Taylor1; S: in Scalar; P2: in out Taylor1) is
  begin
    SetZero(P2);
    Val(P1,S,P2.C(0));
  end Val;

  procedure ValBall(P: in Taylor1; R: in Radius; S: in out Scalar) is
    --- assumes Flt operations are rounded up
    C: Poly1 renames P.C;
    N: Flt := Zero;
  begin
    if P.F=FConst then
      if Check_Consistency and then NotConst(C) then raise Inconsistent_Data; end if;
      Copy(C(0),S);
    elsif (R <= P.R) or else (P.F>PDeg) then
      for I in reverse 1 .. PDeg loop
        N := R*(N+MaxNorm(C(I)));
      end loop;
      BallAt0(N,S);
      Add(P.C(0),S);
    else
      raise Undefined with "Taylor1s.ValBall error: argument too large";
    end if;
  end ValBall;

  procedure Ran(P: in Taylor1; S: in out Scalar) is
  begin
    ValBall(P,P.R,S);
  end Ran;

  function NotSharp(C: Poly1; F,L: Integer) return Boolean is
  begin
    if Not_STrunc then
      for I in F .. L loop
        if not IsSharp(C(I)) then return True; end if;
      end loop;
    end if;
    return False;
  end NotSharp;

  function MaxNorm(R: Radius; P: Taylor1) return Radius is
    --- assumes Flt operations are rounded up
    C: Poly1 renames P.C;
    N: Flt := Zero;
  begin
    if (P.F <= PDeg) and then R>P.R then
      Show2("P.R R",P.R,R);
      raise Domain_Error;
    end if;
    for I in reverse 0 .. EffDeg(C) loop
      N := R*N+MaxNorm(C(I));
    end loop;
    return N;
  end MaxNorm;

  function MaxDist(R: Radius; P1,P2: Taylor1) return Radius is
    --- assumes Flt operations are rounded up
    C1: Poly1 renames P1.C;
    C2: Poly1 renames P2.C;
    N: Flt := Zero;
    S: Scalar;
  begin
    if (P1.F <= PDeg) and then R>P1.R then raise Domain_Error; end if;
    if (P2.F <= PDeg) and then R>P2.R then raise Domain_Error; end if;
    for I in reverse 0 .. IMax(EffDeg(C1),EffDeg(C2)) loop
      Diff(C1(I),C2(I),S);
      N := R*N+MaxNorm(S);
    end loop;
    return N;
  end MaxDist;

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

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

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

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

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

  procedure Compose(C: in SP.Polynom1; P1: in Taylor1; P2: in out Taylor1) is
    D: constant Integer := EffDeg(C);
  begin
    SetZero(P2.C);
    P2.F := P1.F;
    P2.R := P1.R;
    Add(C(D),P2.C(0));
    for I in reverse 0 .. D-1 loop
      Mult(P1,P2);
      Add(C(I),P2.C(0));
    end loop;
  end Compose;

  procedure Reflect(P: in out Taylor1; Offset: in Parity) is
    C: Poly1 renames P.C;
    D: constant Integer := EffDeg(C);
    I: Integer := Offset;
  begin
    while I <= D loop
      Neg(C(I));
      I := I+2;
    end loop;
  end Reflect;

  procedure CutParity(P: in out Taylor1; Offset: in Parity) is
    InvR: constant Radius := One/P.R;
    C: Poly1 renames P.C;
    I: Integer := Offset;
    D: Integer;
  begin
    AdjustF(P,D);
    while I <= D loop
      ModCenter(C(I));
      if I >= P.F and then I<PDeg then
        AddProd(InvR,C(I),C(I+1));
        SetZero(C(I));
        if I=P.F then P.F := P.F+1; end if;
      end if;
      I := I+2;
    end loop;
  end CutParity;

  function ParityErr(P: Taylor1; Offset: Parity) return Radius is
    R: Radius renames P.R;
    C:  Poly1 renames P.C;
    D: constant Integer := EffDeg(C);
    Par: Boolean := (D mod 2)=Offset;
    E: Radius := Zero;
  begin
    for I in reverse 0 .. D loop
      if Par then
        E := R*E+MaxNorm(C(I));
        Par := False;
      else
        E := R*E;
        Par := True;
      end if;
    end loop;
    return E;
  end ParityErr;

  procedure Scale(S: in Scalar; P: in out Taylor1) is
    C:  Poly1 renames P.C;
    R: constant Radius := MaxNorm(S);
  begin
    if R=Zero then
      for I in 1 .. PDeg loop SetZero(C(I)); end loop;
      P.F := FConst;
    else
      P.R := P.R/R;
      if P.F /= FConst then Scale(S,P.C); end if;
    end if;
  end Scale;

  procedure Scale(P1: in Taylor1; S: in Scalar; P2: in out Taylor1) is
    R: constant Radius := MaxNorm(S);
  begin
    if R=Zero then
      SetZero(P2);
      Copy(P1.C(0),P2.C(0));
    else
      if P1.F=FConst then
        Copy(P1,P2);
      else
        Scale(P1.C,S,P2.C);
        P2.F := P1.F;
      end if;
      P2.R := P1.R/R;
    end if;
  end Scale;

  procedure Translate(P1: in Taylor1; S: in Scalar; P2: in out Taylor1) is
    --- assumes Flt operations are rounded up
    A: constant Radius := MaxNorm(S);
    Tmp: LT_Pointer;
  begin
    if A=Zero then
      Copy(P1,P2);
    elsif P1.F=FConst then
      if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
      Copy(P1,P2);
    else
      TPool.Allocate(Tmp);
      Translate(P1.C,S,P2.C,Tmp.all.Data.C);
      TPool.Recycle(Tmp);
      if P1.F>PDeg then
        P2.F := PDeg1;
        P2.R := P1.R;
      else
        P2.F := 0;
        P2.R := -(A-P1.R);
        LiftErrs(P2,PDeg);
      end if;
    end if;
  end Translate;

  procedure Translate(P1: in Taylor1; S: in Scalar; P2: in out Taylor1; R: in Radius) is
  begin
    Translate(P1,S,P2);
    SetRho(R,P2);
  end Translate;

  procedure Der(R: in Radius; P: in out Taylor1) is
    C: Poly1 renames P.C;
  begin
    if P.F>PDeg then
      if P.F=FConst then
        if Check_Consistency and then NotConst(C) then raise Inconsistent_Data; end if;
        SetZero(C);
      else
        Num_Der(C);
        P.F := TrueF(P);
      end if;
      P.R := R;
      return;
    end if;
    if P.R <= R then
      Show2("",P.R,R);
      raise Domain_Error with Show0("Taylors1.Der: need P.R>R");
    end if;
    LiftErrs(P,1);
    declare
      S: constant Scalar := Inv(Log(Scal(P.R)/R));
      E1: constant Scalar := Exp(Scal(-1));
      RM: constant Flt := MaxNorm(S);
      FM: constant Flt := MaxNorm(S*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;
    begin
      for I in P.F .. PDeg loop
        ErrMult(DerFac(I),C(I));
      end loop;
    end;
    Num_Der(C);
    P.F := P.F-1;
    P.R := R;
  end Der;

  procedure Num_ValDer(P: in Poly1; S: in Scalar; C: in out Scalar) is
    D: constant Natural := EffLast(P);
    Tmp: Scalar;
  begin
    Prod(Flt(D),P(D),C);
    for N in reverse 1 .. D-1 loop
      Mult(S,C);
      AddProd(Flt(N),P(N),C,Tmp);
    end loop;
  end Num_ValDer;

  procedure ValDer(P: in Taylor1; S: in Scalar; D: in out Scalar) is
    Tmp: LT_Pointer;
  begin
    if P.F>PDeg then
      if P.F=FConst then
        if Check_Consistency and then NotConst(P.C) then raise Inconsistent_Data; end if;
        SetZero(D);
      else
        Num_ValDer(P.C,S,D);
      end if;
    else
      TPool.Allocate(Tmp);
      declare
        R: constant Radius := MaxNorm(S);
        Q: Taylor1 renames Tmp.all.Data;
      begin
        Copy(P,Q);
        if R=Zero then
          LiftErrs(Q,1);
          Copy(Q.C(1),D);
        else
          Der(R,Q);
          Val(Q,S,D);
        end if;
      end;
      TPool.Recycle(Tmp);
    end if;
  end ValDer;

  procedure DivArg(P: in out Taylor1) is
    --- [P(x)-P(0)]/x
    C: Poly1 renames P.C;
  begin
    if P.F=FConst then
      if Check_Consistency and then NotConst(C) then raise Inconsistent_Data; end if;
      SetZero(C(0));
    else
      LiftErrs(P,1);
      for I in 1 .. PDeg loop Copy(C(I),C(I-1)); end loop;
      SetZero(C(PDeg));
      if P.F <= PDeg then P.F := P.F-1; end if;
    end if;
  end DivArg;

  procedure AntiDer(P: in out Taylor1) is
    C: Poly1 renames P.C;
    E: Scalar;
  begin
    ToErr(C(PDeg),E);
    pragma Warnings (Off);
    for I in reverse 2 .. EffDeg(C) loop
      Quot(C(I-1),Flt(I),C(I));
    end loop;
    pragma Warnings (On);
    Copy(C(0),C(1));
    if not IsZero(E) then
      Div(Flt(PDeg+1),E);
      AddProd(P.R,E,C(PDeg));
    end if;
    SetZero(C(0));
    if P.F<PDeg then P.F := P.F+1; end if;
  end AntiDer;

  procedure AntiDer(P1: in Taylor1; P2: in out Taylor1) is
    C1: Poly1 renames P1.C;
    C2: Poly1 renames P2.C;
    E: Scalar renames C2(0);
  begin
    for I in reverse 2 .. PDeg loop Quot(C1(I-1),Flt(I),C2(I)); end loop;
    Copy(C1(0),C2(1));
    if Not_STrunc then
      ToErr(C1(PDeg),E);
      Div(Flt(PDeg+1),E);
      Add(E,C2(PDeg));
    end if;
    SetZero(C2(0));
    if P1.F<PDeg then P2.F := P1.F+1; else P2.F := P1.F; end if;
    P2.R := P1.R;
  end AntiDer;

  function Integral(P: Taylor1; S1,S2: Scalar) return Scalar is
    Tmp: LT_Pointer;
    V1,V2: Scalar;
  begin
    TPool.Allocate(Tmp);
    AntiDer(P,Tmp.all.Data);
    Val(Tmp.all.Data,S1,V1);
    Val(Tmp.all.Data,S2,V2);
    TPool.Recycle(Tmp);
    Sub(V1,V2);
    return V2;
  end Integral;

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

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

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

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

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

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

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

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

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

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

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

  procedure Inv(P1: in Taylor1; P2: in out Taylor1) is
  begin
    if P1.F=FConst then
      if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
      SetZero(P2);
      Inv(P1.C(0),P2.C(0));
    else
      declare
        Iter: constant Integer := Choose(STrunc,PDeg,IMax(PDeg+8,32));
      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 Taylor1; P2: in out Taylor1) is
  begin
    if P1.F=FConst then
      if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
      SetZero(P2);
      Sqrt(P1.C(0),P2.C(0));
    else
      declare
        Iter: constant Integer := Choose(STrunc,PDeg,IMax(PDeg+8,32));
      begin
        TFS.Series_QPower(IRat(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 Taylor1; P2: in out Taylor1) is
  begin
    if K=1 then
      Copy(P1,P2);
    elsif P1.F=FConst then
      if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
      SetZero(P2);
      Root(K,P1.C(0),P2.C(0));
    else
      declare
        Iter: constant Integer := Choose(STrunc,PDeg,IMax(PDeg+8,32));
      begin
        TFS.Series_QPower(IRat(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 Taylor1; P2: in out Taylor1) is
  begin
    if P1.F=FConst then
      if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
      SetZero(P2);
      QPower(Q,P1.C(0),P2.C(0));
    elsif IDen(Q)=1 then
      IPower(INum(Q),P1,P2);
    else
      declare
        Iter: constant Integer := Choose(STrunc,PDeg,IMax(PDeg+8,32));
      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 Taylor1; P2: in out Taylor1) is
  begin
    if P1.F=FConst then
      if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
      SetZero(P2);
      Exp(P1.C(0),P2.C(0));
    else
      declare
        Iter: constant Integer := Choose(STrunc,PDeg,IMax(PDeg+8,32));
      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 Taylor1; PC,PS: in out Taylor1) is
  begin
    if P.F=FConst then
      if Check_Consistency and then NotConst(P.C) then raise Inconsistent_Data; end if;
      SetZero(PC);
      Cos(P.C(0),PC.C(0));
      SetZero(PS);
      Sin(P.C(0),PS.C(0));
    else
      declare
        Iter: constant Integer := Choose(STrunc,PDeg/2,IMax(PDeg+8,32)/2); -- 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 Taylor1; P2: in out Taylor1) is
  begin
    if P1.F=FConst then
      if Check_Consistency and then NotConst(P1.C) then raise Inconsistent_Data; end if;
      SetZero(P2);
      Log(P1.C(0),P2.C(0));
    else
      declare
        Iter: constant Integer := Choose(STrunc,PDeg,IMax(PDeg+8,32));
      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 Taylor1; P2: in out Taylor1) is
    --- approx and modulo a constant
    Tmp: LT_Pointer;
  begin
    TPool.Allocate(Tmp);
    declare
      Q: Taylor1 renames Tmp.all.Data;
    begin
      Prod(P1,P1,Q);
      Center(Q);
      Neg(Q);
      Add(1,Q);
      QPower(IRat(-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;
    TPool.Recycle(Tmp);
    Center(P2);
    Inv_MDer(P2.C,False);   -- integrate P2
  end ArcSin_Part;

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

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

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

  function ExpTail(K: Natural; R: Radius) return Radius is
    --- bound on R/[(K+1)!]+R^2/[(K+2)!]+R^3/[(K+3)!]+...
    Q: Flt := R/Flt(K+1);
    E: Flt;
  begin
    if Q<Half then
      E := Q/abs(Q-One);
    else
      E := Q;
      declare
        R0: constant Positive := 4096;
        EN: Flt := Q;
      begin
        if R>Flt(R0) then raise Sorry; end if;
        for N in K+2 .. 2*R0 loop
          Q := R/Flt(N);
          EN := Q*EN;
          exit when Q<Half;
          E := E+EN;
        end loop;
        E := E+EN/abs(Q-One);
      end;
    end if;
    for N in 2 .. K loop E := E/Flt(N); end loop;
    return E;
  end ExpTail;

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

  procedure CosTaylor(S: in Scalar; R: in Radius; P: in out Taylor1) is
    C: Poly1 renames P.C;
  begin
    Cos(S,C(0));
    Sin(S,C(1));
    Neg(C(1));
    for N in 2 .. PDeg loop
      Quot(C(N-2),Flt(N*(1-N)),C(N));
    end loop;
    P.R := R;
    if STrunc then
      P.F := PDeg1;
    else
      if not IsReal(S) then
        raise Not_Implemented; -- bound below is only for real S
      end if;
      P.F := PDeg;
      Add(BallAt0(ExpTail(PDeg,R)),C(PDeg));
    end if;
  end CosTaylor;

  procedure CosTaylor(K: in Integer; S: in Scalar; R: in Radius; P: in out Taylor1) is
  begin
    CosTaylor(Scal(K),S,R,P);
  end CosTaylor;

  procedure CosTaylor(K,S: in Scalar; R: in Radius; P: in out Taylor1) is
    C: Poly1 renames P.C;
    RK: constant Flt := Sup(abs(K));
    T: Scalar;
  begin
    if IsZero(K) then
      SetOne(R,P);
    else
      Prod(K,S,T);
      CosTaylor(T,RK*R,P);
      Copy(SOne,T);
      for N in 1 .. PDeg loop
        Mult(K,T);
        Mult(T,C(N));
      end loop;
      P.R := R;
    end if;
  end CosTaylor;

  procedure SinTaylor(S: in Scalar; R: in Radius; P: in out Taylor1) is
    C: Poly1 renames P.C;
  begin
    Sin(S,C(0));
    Cos(S,C(1));
    for N in 2 .. PDeg loop
      Quot(C(N-2),Flt(N*(1-N)),C(N));
    end loop;
    P.R := R;
    if STrunc then
      P.F := PDeg1;
    else
      if not IsReal(S) then
        raise Not_Implemented; -- bound below is only for real S
      end if;
      P.F := PDeg;
      Add(BallAt0(ExpTail(PDeg,R)),C(PDeg));
    end if;
  end SinTaylor;

  procedure SinTaylor(K: in Integer; S: in Scalar; R: in Radius; P: in out Taylor1) is
  begin
    SinTaylor(Scal(K),S,R,P);
  end SinTaylor;

  procedure SinTaylor(K,S: in Scalar; R: in Radius; P: in out Taylor1) is
    C: Poly1 renames P.C;
    RK: constant Flt := Sup(abs(K));
    T: Scalar;
  begin
    if IsZero(K) then
      SetZero(R,P);
    else
      Prod(K,S,T);
      SinTaylor(T,RK*R,P);
      Copy(SOne,T);
      for N in 1 .. PDeg loop
        Mult(K,T);
        Mult(T,C(N));
      end loop;
      P.R := R;
    end if;
  end SinTaylor;

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

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

  procedure Roots3(B,C,D: in Taylor1; U0,U1,U2,V: in out Taylor1) 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),C.C(0),D.C(0),U0.C(0),U1.C(0),U2.C(0),V.C(0));
    else
      TR.NewtonRoots3(B,C,D,U0,U1,U2,V,Steps);
    end if;
  end Roots3;

  ------------------ FBalls

  function MaxNorm(B: in FBall) return Radius is
  begin
    return abs(B.C)+B.R;
  end MaxNorm;

  procedure Enclose(R1,R2: in Flt; B: out FBall) is
    H1:   constant Flt := Half*R1;
    H2:   constant Flt := Half*R2;
    NegC: constant Flt := -H1-H2;
  begin
    B.C := H1+H2;
    if H1 <= H2 then
      B.R := (H2-H1)+(B.C+NegC);
    else
      B.R := (H1-H2)+(B.C+NegC);
    end if;
  end Enclose;

  procedure AddProd(B1,B2: in FBall; B3: in out FBall) is
    C: Flt renames B3.C;
    NegC: constant Flt := (-B1.C)*B2.C-C;
  begin
    C := B1.C*B2.C+C;
    B3.R := abs(B1.C)*B2.R+B1.R*(abs(B2.C)+B2.R)+B3.R+(C+NegC);
  end AddProd;

  ---------------------- numerical

  function GuessRho(P: Taylor1) return Radius is
    use Flt_EF;
    DMin: constant Integer := PDeg/4;
    DMax: constant Integer := (3*PDeg)/4;
    C: Poly1 renames P.C;
    I: Integer := DMin;
    N: Integer := 0;
    LogQ: Flt := Zero;
    Num: Flt;
    pragma Warnings (Off);
    Den: Flt := MaxNorm(C(I-2))+MaxNorm(C(I-1));
    pragma Warnings (On);
  begin
    while I<DMax loop
      Num := Den;
      Den := MaxNorm(C(I))+MaxNorm(C(I+1));
      if not (Num=Zero or else Den=Zero) then
        LogQ := LogQ+Log(Num/Den);
        N := N+2;
      end if;
      I := I+2;
    end loop;
    return Exp(LogQ/Flt(N));
  exception
    when CONSTRAINT_ERROR =>
      Show0("GuessRho: caught CONSTRAINT_ERROR");
      return Zero;
  end GuessRho;

  function GuessRhoInv(P: Taylor1) return Radius is
    Tmp: LT_Pointer;
    R: Radius := Zero;
  begin
    TPool.Allocate(Tmp);
    declare
      Q: Taylor1 renames Tmp.all.Data;
    begin
      Inv(P,Q);
      R := GuessRho(Q);
    exception
      when others => null;
    end;
    TPool.Recycle(Tmp);
    return R;
  end GuessRhoInv;

  procedure NumInterpolate(P0: in Taylor1; P: in out Taylor1; Fac: in Flt) is
    --- P := (1-Fac)*P0+Fac*P
  begin
    if Fac /= One then
      if Fac /= Zero then
        declare
          Fac0: constant Flt := One-Fac;
          C0: Poly1 renames P0.C;
          C: Poly1 renames P.C;
          Tmp: Scalar;
        begin
          for I in Power loop
            Mult(Fac,C(I),Tmp);
            AddProd(Fac0,C0(I),C(I),Tmp);
          end loop;
        end;
      else
        Copy(P0,P);
      end if;
    end if;
  end NumInterpolate;

  procedure NumTrunc(P: in out Taylor1; Deg: in Power := PDeg; Cut: in Radius := Zero) is
    PosCut: constant Boolean := Cut>Zero;
    C: Poly1 renames P.C;
    procedure Truncate(S: in out Scalar) is
    begin
      Center(S);
      if PosCut and then MaxNorm(S)<Cut then SetZero(S); end if;
    end Truncate;
  begin
    Truncate(C(0));
    if P.F /= FConst then
      for I in 1 .. Deg loop Truncate(C(I)); end loop;
      for I in Deg+1 .. PDeg loop SetZero(C(I)); end loop;
      P.F := PDeg1;
    end if;
  end NumTrunc;

  procedure NumIProd(P1,P2: in Taylor1; S: in out Scalar) is
    C1: Poly1 renames P1.C;
    C2: Poly1 renames P2.C;
  begin
    SetZero(S);
    for I in reverse Power loop AddProd(C1(I),C2(I),S); end loop;
  end NumIProd;

begin

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

end Taylors1;
