File : taylors-ops.adb

with Messages;
use Messages;

pragma Elaborate_All (Messages);

package body Taylors.Ops is

  ScalOne:   constant Scalar := Scal(One);

  function IMin(I1,I2: Integer) return Integer renames Integer'Min;
  function IMax(I1,I2: Integer) return Integer renames Integer'Max;

  pragma Inline (Imin,Imax);

  procedure ComponentSplit(K: in Power; P: in Taylor; P0: out Taylor; E: out Scalar) is
  begin
    if K>P.D then
      P0.D := -1;
      E := ScalZero;
      return;
    end if;
    P0.D := 0;
    if K<Dho then
      P0.C(0) := P.C(K);
      E := ScalZero;
    else
      Split(P.C(K),P0.C(0),E);
    end if;
    if Trunc then return; end if;
    P0.E := ScalZero;
    if K=0 then E := E+P.E; end if;
  end ComponentSplit;

  procedure AddComponent(D: in Power; S: in Scalar; P: in out Taylor) is
  begin
    if S=ScalZero then return; end if;
    if D > P.D then
      for K in P.D+1 .. D-1 loop
        P.C(K) := ScalZero;
      end loop;
      P.C(D) := S;
      if P.D<0 then P.E := ScalZero; end if;
      P.D := D;
    else
      P.C(D) := P.C(D)+S;
    end if;
  end AddComponent;

  procedure MultComponent(D: in Power; S: in Scalar; P: in out Taylor) is
  begin
    if D <= P.D then
      P.C(D) := S*P.C(D);
    end if;
  end MultComponent;

  procedure AddErrComp(S: in Scalar; P: in out Taylor) is
  begin
    if Trunc or else S=ScalZero then return; end if;
    if P.D >= 0 then
      P.E := P.E+Abs(S);
    else
      P.C(0) := ScalZero;
      P.D := 0;
      P.E := Abs(S);
    end if;
  end AddErrComp;

  procedure MultErrComp(S: in Scalar; P: in out Taylor) is
  begin
    if Trunc or else P.D<0 then return; end if;
    P.E := Abs(S)*P.E;
  end MultErrComp;

  function Coeff(R: Radius; K: Power; P: Taylor) return Scalar is
  begin
    if P.D<0 then
      return ScalZero;
    elsif Trunc then
      return Component(K,P);
    elsif K <= Dho then
      return Component(K,P)+Ball(P.E)/Scal(R)**K;
    end if;
    declare
      K1: constant Integer := IMin(K-1,P.D);
      E: Scalar := Ball(P.E)/Scal(R)**(Dho-1);
    begin
      for N in Dho .. K1 loop
        E := E/R+ResetCenter(P.C(N));
      end loop;
      return Component(K,P)+E/Scal(R)**(K-K1);
    end;
  end Coeff;

  procedure ResetCenter(P: in out Taylor; Kmin: in Natural := 0; Kmax: in Integer := Size) is
  begin
    if Trunc then
      ResetComponent(P,Kmin,Kmax);
    else
      for K in Kmin .. Imin(Kmax,P.D) loop
        ResetCenter(P.C(K));
      end loop;
    end if;
  end ResetCenter;

  function Val(R: Radius; Z: Scalar; P: Taylor) return Scalar is
    S: Scalar := ScalZero;
  begin
    for K in reverse 0 .. P.D loop
      S := Z*S+P.C(K);
    end loop;
    if Trunc or else P.D<0 then return S; end if;
    if SupAbs(Z)>R then
      Message("Taylors.Ops.Val",Domain_Violation);
    end if;
    return S+Ball(P.E);
  end Val;

  function Norm(R: Radius; P: Taylor) return Scalar is
    E: Scalar := ScalZero;
  begin
    for K in reverse 0 .. P.D loop
      E := R*E+Abs(P.C(K));
    end loop;
    if Trunc or else P.D<0 then return E; end if;
    return E+P.E;
  end Norm;

  procedure ErrNorm(R: in Radius; Pw,P1: in Taylor; P2: out Taylor; E: out Scalar) is
--- splits P1 into P2+P, E is weighted L1 norm of P
    D: constant Degree := Imin(Pw.D,P1.D);
    F: Scalar := ScalZero;
  begin
    E := ScalZero;
    if D<0 or else Trunc then
      Copy(P1,P2);
      return;
    end if;

    for K in reverse D+1 .. P1.D loop
      P2.C(K) := P1.C(K);
    end loop;

    for K in reverse D+1 .. Pw.D loop
      F := Max(F,Pw.C(K));
    end loop;

    for K in reverse Dho .. D loop
      P2.C(K) := Center(P1.C(K));
      F := Max(F,Pw.C(K));
      E := R*E+Abs(F*ResetCenter(P1.C(K)));   -- add to E only when weight is nonzero
    end loop;

    for K in reverse 0 .. Imin(Dho-1,D) loop
      P2.C(K) := P1.C(K);
      F := Max(F,Pw.C(K));
      E := R*E;
    end loop;

    P2.D := P1.D;
    P2.E := ScalZero;
    E := E+F*P1.E;
  end ErrNorm;

  procedure Neg(P1: in Taylor; P2: out Taylor) is
  begin
    P2.D := P1.D;
    for K in 0 .. P1.D loop
      P2.C(K) := -P1.C(K);
    end loop;
    P2.E := P1.E;                  -- garbage if Trunc
  end Neg;

  function "-"(P: Taylor) return Taylor is
    Pt: Taylor;
  begin
    Neg(P,Pt);
    return Pt;
  end "-";

  procedure Mult(S: in Scalar; P: in out Taylor) is
  begin
    for K in 0 .. P.D loop
      P.C(K) := S*P.C(K);
    end loop;
    if Trunc or else P.D<0 then return; end if;
    P.E := Abs(S)*P.E;
  end Mult;

  procedure Prod(S: in Scalar; P1: in Taylor; P2: out Taylor) is
  begin
    P2.D := P1.D;
    for K in 0 .. P1.D loop
      P2.C(K) := S*P1.C(K);
    end loop;
    if Trunc or else P1.D<0 then return; end if;
    P2.E := Abs(S)*P1.E;
  end Prod;

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

  procedure Add(P1: in Taylor; P2: in out Taylor) is
  begin
    if P1.D <= P2.D then
      if P1.D<0 then return; end if;
      for K in 0 .. P1.D loop
        P2.C(K) := P2.C(K)+P1.C(K);
      end loop;
    else
      if P2.D<0 then Copy(P1,P2); return; end if;
      for K in 0 .. P2.D loop
        P2.C(K) := P2.C(K)+P1.C(K);
      end loop;
      for K in P2.D+1 .. P1.D loop
        P2.C(K) := P1.C(K);
      end loop;
      P2.D := P1.D;
    end if;
    if Trunc then return; end if;
    P2.E := P2.E+P1.E;
  end Add;

  procedure Sum(P1,P2: in Taylor; P3: out Taylor) is
  begin
    if P1.D <= P2.D then
      if P1.D<0 then Copy(P2,P3); return; end if;
      for K in 0 .. P1.D loop
        P3.C(K) := P1.C(K)+P2.C(K);
      end loop;
      for K in P1.D+1 .. P2.D loop
        P3.C(K) := P2.C(K);
      end loop;
      P3.D := P2.D;
    else
      if P2.D<0 then Copy(P1,P3); return; end if;
      for K in 0 .. P2.D loop
        P3.C(K) := P1.C(K)+P2.C(K);
      end loop;
      for K in P2.D+1 .. P1.D loop
        P3.C(K) := P1.C(K);
      end loop;
      P3.D := P1.D;
    end if;
    if Trunc then return; end if;
    P3.E := P1.E+P2.E;
  end Sum;

  function "+"(P1,P2: Taylor) return Taylor is
    Pt: Taylor;
  begin
    Sum(P1,P2,Pt);
    return Pt;
  end "+";

  procedure Sub(P1: in Taylor; P2: in out Taylor) is
  begin
    if P1.D <= P2.D then
      if P1.D<0 then return; end if;
      for K in 0 .. P1.D loop
        P2.C(K) := P2.C(K)-P1.C(K);
      end loop;
    else
      if P2.D<0 then Neg(P1,P2); return; end if;
      for K in 0 .. P2.D loop
        P2.C(K) := P2.C(K)-P1.C(K);
      end loop;
      for K in P2.D+1 .. P1.D loop
        P2.C(K) := -P1.C(K);
      end loop;
      P2.D := P1.D;
    end if;
    if Trunc then return; end if;
    P2.E := P2.E+P1.E;
  end Sub;

  procedure Diff(P1,P2: in Taylor; P3: out Taylor) is
  begin
    if P1.D <= P2.D then
      if P1.D<0 then Neg(P2,P3); return; end if;
      for K in 0 .. P1.D loop
        P3.C(K) := P1.C(K)-P2.C(K);
      end loop;
      for K in P1.D+1 .. P2.D loop
        P3.C(K) := -P2.C(K);
      end loop;
      P3.D := P2.D;
    else
      if P2.D<0 then Copy(P1,P3); return; end if;
      for K in 0 .. P2.D loop
        P3.C(K) := P1.C(K)-P2.C(K);
      end loop;
      for K in P2.D+1 .. P1.D loop
        P3.C(K) := P1.C(K);
      end loop;
      P3.D := P1.D;
    end if;
    if Trunc then return; end if;
    P3.E := P1.E+P2.E;
  end Diff;

  function "-"(P1,P2: Taylor) return Taylor is
    Pt: Taylor;
  begin
    Diff(P1,P2,Pt);
    return Pt;
  end "-";

  procedure MultAdd(S: in Scalar; P1: in Taylor; P2: in out Taylor) is
  begin
    if P1.D <= P2.D then
      if P1.D<0 or else S=ScalZero then return; end if;
      for K in 0 .. P1.D loop
        P2.C(K) := P2.C(K)+S*P1.C(K);
      end loop;
    else
      if P2.D<0 then Prod(S,P1,P2); return; end if;
      if S=ScalZero then return; end if;
      for K in 0 .. P2.D loop
        P2.C(K) := P2.C(K)+S*P1.C(K);
      end loop;
      for K in P2.D+1 .. P1.D loop
        P2.C(K) := S*P1.C(K);
      end loop;
      P2.D := P1.D;
    end if;
    if Trunc then return; end if;
    P2.E := P2.E+Abs(S)*P1.E;
  end Multadd;

  procedure Scale(S: in Scalar; P1: in Taylor; P2: out Taylor) is
  begin
    P2.D := P1.D;
    if P1.D<0 then return; end if;
    P2.C(0) := P1.C(0)/S;
    declare
      Sp: Scalar := ScalOne;
    begin
      for I in 1 .. P2.D loop
        P2.C(I) := P1.C(I)*Sp;
        Sp := Sp*S;
      end loop;
    end;
    if Trunc then return; end if;
    P2.E := P1.E/Abs(S);
  end Scale;

  procedure Scale(S: in Scalar; P: in out Taylor) is
  begin
    if P.D<0 then return; end if;
    P.C(0) := P.C(0)/S;
    declare
      Sp: Scalar := S;
    begin
      for I in 2 .. P.D loop
        P.C(I) := P.C(I)*Sp;
        Sp := Sp*S;
      end loop;
    end;
    if Trunc then return; end if;
    P.E := P.E/Abs(S);
  end Scale;

  function Scale(S: Scalar; P: Taylor) return Taylor is
    Pt: Taylor;
  begin
    Scale(S,P,Pt);
    return Pt;
  end Scale;

  procedure ScaleRem1(Eps: in Scalar; P1: in Taylor; P2: out Taylor; IsPoly: in Boolean := False) is
  begin
    if P1.D<0 then P2.D := -1; return; end if;
    if not IsPoly then
      Message("Taylors.Ops.ScaleRem1",Not_Implemented);
    end if;
    P2.E := ScalZero;
    P2.C(0) := Eps*P1.C(0);
    if P1.D <= 1 then P2.D := 0; return; end if;
    P2.C(1) := ScalZero;
    declare
      S: constant Scalar := Inv(ScalOne+Eps);
      Fac: Scalar := ScalZero;
    begin
      for K in 2 .. P1.D loop
        Fac := (Fac-Eps)*S;
        P2.C(K) := Fac*P1.C(K);
      end loop;
    end;
    P2.D := P1.D;
  end ScaleRem1;

  procedure ScaleRem2(Eps: in Scalar; P1: in Taylor; P2: out Taylor; IsPoly: in Boolean := False) is
  begin
    if P1.D <= 1 then P2.D := -1; return; end if;
    if not IsPoly then
      Message("Taylors.Ops.ScaleRem2",Not_Implemented);
    end if;
    P2.E := ScalZero;
    P2.C(0) := ScalZero;
    P2.C(1) := ScalZero;
    declare
      Eps2: constant Scalar := Sqr(Eps);
      S: constant Scalar := Inv(ScalOne+Eps);
      Fac: Scalar := ScalZero;
    begin
      for K in 2 .. P1.D loop
        Fac := (Fac+Rep(K-1)*Eps2)*S;
        P2.C(K) := Fac*P1.C(K);
      end loop;
    end;
    P2.D := P1.D;
  end ScaleRem2;

  procedure Fmax1(R: in Scalar; M: out Natural; F: out Scalar) is
    X: constant Scalar := Inv(-Log(R));
  begin
    M := Integer(Rep'Ceiling(Sup(X)));
    F := X*Exp(-ScalOne);
  end Fmax1;

  procedure Der(R1,R2: in Radius; P1: in Taylor; P2: out Taylor; IsPoly: in Boolean := False) is
--  error in P1.C(Dho) effectively doubles, but degree info is kept
  begin
    if P1.D<0 then P2.D := -1; return; end if;
    for K in 1 .. P1.D loop
      P2.C(K-1) := Rep(K)*P1.C(K);
    end loop;
    P2.D := P1.D-1;
    if Trunc then return; end if;
    P2.E := ScalZero;
    if IsPoly then return; end if;
    if R2>=R1 then
      Message("Taylors.Ops.Der",Domain_Violation);
    end if;
    declare
      Ri: constant Scalar := Scal(R1)/R2;
      M: Natural;
      E,F: Scalar;
    begin
      Fmax1(Scal(R2)/R1,M,F);
      P2.E := P1.E*F/R2;
      if P1.D <= Dho then
        if P1.D<Dho then return; end if;
        P2.D := Dho;
        P2.C(Dho) := ScalZero;
      end if;
      F := F*Ri**Dho;
      E := ResetCenter(P1.C(Dho))*F/R2;
      for K in Dho+1 .. IMin(M,P1.D) loop
        F := F*Ri;
        ErrMult(SupAbs(F/Rep(K)),P2.C(K-1));
      end loop;
      P2.C(Dho) := P2.C(Dho)+E;
    end;
  end Der;

  function Der(R1,R2: Radius; P: Taylor; IsPoly: in Boolean := False) return Taylor is
    Pt: Taylor;
  begin
    Der(R1,R2,P,Pt,IsPoly);
    return Pt;
  end Der;

  function DerVal(R: Radius; Z: Scalar; P: Taylor) return Scalar is
    Pd: Taylor;
  begin
    Der(R,SupAbs(Z),P,Pd);
    return Val(SupAbs(Z),Z,Pd);
  end DerVal;

  procedure DerMult(R,Fmax: in Radius; S: in Scalar; P: in out Taylor; Plain: in Boolean := False) is
--  error in P.C(Dho) effectively doubles, but degree info is kept
  begin
    if P.D<0 then return; end if;
    if Plain then
      for K in 1 .. P.D loop
        P.C(K-1) := Rep(K)*S*P.C(K);
      end loop;
      P.D := P.D-1;
      return;
    end if;

    P.E := Fmax*P.E/R;                   -- operator_norm <= Fmax/R
    if P.D<Dho then
      if P.D=0 then
        P.C(0) := ScalZero;
        return;
      end if;
      for K in 1 .. P.D loop
        P.C(K-1) := Cap(Fmax,Rep(K)*S)*P.C(K);
      end loop;
      P.D := P.D-1;
      return;
    end if;

    declare
      E: constant Scalar := Fmax*ResetCenter(P.C(Dho))/R;
    begin
      for K in 1 .. Dho loop
        P.C(K-1) := Cap(Fmax,Rep(K)*S)*P.C(K);
      end loop;
      if P.D=Dho then
        P.C(Dho) := E;
        return;
      end if;
      for K in Dho+1 .. P.D loop
        P.C(K-1) := Cap(Fmax,Rep(K)*S)*Center(P.C(K))+Fmax*ResetCenter(P.C(K));
      end loop;
      P.C(Dho) := P.C(Dho)+E;
      P.D := P.D-1;
    end;
  end DerMult;

  procedure DerProd(R,Fmax: in Radius; S: in Scalar; P1: in Taylor; P2: out Taylor; Plain: in Boolean := False) is
--  error in P.C(Dho) effectively doubles, but degree info is kept
  begin
    P2.D := P1.D;
    if P2.D<0 then return; end if;
    if Plain then
      for K in 1 .. P2.D loop
        P2.C(K-1) := Rep(K)*S*P1.C(K);
      end loop;
      P2.D := P2.D-1;
      P2.E := ScalZero;
      return;
    end if;

    P2.E := Fmax*P1.E/R;                 -- operator_norm <= Fmax/R
    if P2.D<Dho then
      if P2.D=0 then
        P2.C(0) := ScalZero;
        return;
      end if;
      for K in 1 .. P2.D loop
        P2.C(K-1) := Cap(Fmax,Rep(K)*S)*P1.C(K);
      end loop;
      P2.D := P2.D-1;
      return;
    end if;

    for K in 1 .. Dho loop
      P2.C(K-1) := Cap(Fmax,Rep(K)*S)*P1.C(K);
    end loop;
    if P2.D=Dho then
      P2.C(Dho) := Fmax*ResetCenter(P1.C(Dho))/R;
      return;
    end if;
    for K in Dho+1 .. P2.D loop
      P2.C(K-1) := Cap(Fmax,Rep(K)*S)*Center(P1.C(K))+Fmax*ResetCenter(P1.C(K));
    end loop;
    P2.C(Dho) := P2.C(Dho)+Fmax*ResetCenter(P1.C(Dho))/R;
    P2.D := P2.D-1;
  end DerProd;

  procedure ZxDer(P1: in Taylor; P2: out Taylor; IsPoly: in Boolean := False) is
  begin
    if not IsPoly then
      Message("Taylors.Ops.ZxDer",Not_Implemented);
    end if;
    if P1.D<=0 then
      P2.D := -1;
      return;
    end if;
    P2.D := P1.D;
    P2.C(0) := ScalZero;
    for K in 1 .. P1.D loop
      P2.C(K) := Rep(K)*P1.C(K);
    end loop;
    P2.E := ScalZero;
  end ZxDer;

  procedure DerDer(P: in out Taylor; IsPoly: in Boolean := False) is
  begin
    if not IsPoly then
      Message("Taylors.Ops.DerDer",Not_Implemented);
    end if;
    for K in 2 .. P.D loop
      P.C(K-2) := Rep(K*(K-1))*P.C(K);
    end loop;
    P.D := Imax(P.D-2,-1);
  end DerDer;

  function DerDer(P: Taylor; IsPoly: in Boolean := False) return Taylor is
    Pr: Taylor;
  begin
    if not IsPoly then
      Message("Taylors.Ops.DerDer",Not_Implemented);
    end if;
    Pr.D := Imax(P.D-2,-1);
    for K in 2 .. P.D loop
      Pr.C(K-2) := Rep(K*(K-1))*P.C(K);
    end loop;
    Pr.E := ScalZero;
    return Pr;
  end DerDer;

  procedure DegConvert(R: in Radius; D: in Degree; P: in out Taylor) is
  begin
    if D >= P.D then return; end if;
    if not Trunc then
      declare
        E: Scalar := ScalZero;
      begin
        for K in reverse D+1 .. P.D loop
          E := R*E+Ball(P.C(K));
        end loop;
        if D<Dho then
          if D<0 then
            Message("Taylors.Ops.DegConvert",Domain_Violation);
          end if;
          P.E := P.E+SupAbs(E)*Scal(R)**(D+1);
        else
          P.C(D) := P.C(D)+R*E;
        end if;
      end;
    end if;
    P.D := D;
  end DegConvert;

  procedure XMultAdd(R: in Radius; K: in Power; S: in Scalar; P1: in Taylor; P2: in out Taylor; D: in Degree := Size) is
  begin
    if Trunc then
      if D<0 then P2.D := -1; return; end if;
      if D<K then DegConvert(R,D,P2); return; end if;
    elsif D<Dho then
      XMultAdd(R,K,S,P1,P2,Dho);
      DegConvert(R,D,P2);
      return;
    end if;

    DegConvert(R,D,P2);
    if P1.D<0 or else S=ScalZero then return; end if;

    declare
      D1: constant Natural := K+P1.D;
      D2: Natural := IMin(D1,D);
    begin
      if P2.D >= D2 then                       -- l.o. stuff fits below P2.D
        for I in K .. D2 loop
          P2.C(I) := P2.C(I)+S*P1.C(I-K);
        end loop;
      elsif P2.D >= K then                     -- l.o. stuff above and below P2.D
        for I in K .. P2.D loop
          P2.C(I) := P2.C(I)+S*P1.C(I-K);
        end loop;
        for I in P2.D+1 .. D2 loop
          P2.C(I) := S*P1.C(I-K);
        end loop;
        P2.D := D2;
      else                                     -- no l.o. stuff below P2.D
        for I in P2.D+1 .. IMin(K-1,D2) loop
          P2.C(I) := ScalZero;
        end loop;
        for I in K .. D2 loop
          P2.C(I) := S*P1.C(I-K);
        end loop;
        if P2.D<0 then P2.E := ScalZero; end if;
        P2.D := D2;
      end if;
      if Trunc then return; end if;

      if K<Dho then                            -- add to general error
        P2.E := P2.E+Abs(S)*P1.E*Scal(R)**K;
      else
        D2 := Imin(K,P2.D);
        P2.C(D2) := P2.C(D2)+Abs(S)*P1.E*Scal(R)**(K-D2);
      end if;

      if D1>D then                             -- add h.o. stuff
        D2 := IMax(K,D+1);
        declare
          E: Scalar := ScalZero;
        begin
          for I in reverse D2 .. D1 loop
            E := R*E+Ball(P1.C(I-K));
          end loop;
          P2.C(D) := P2.C(D)+S*E*Scal(R)**(D2-D);
        end;
      end if;
    end;
  end XMultAdd;

  procedure Mult(R: in Radius; P1: in Taylor; P2: in out Taylor; D: in Degree := Size) is
  begin
    if P2.D<0 then return; end if;
    if P1.D<0 then P2.D := -1; return; end if;

    if Trunc then
      declare
        D2: constant Natural := Imin(D,P1.D+P2.D);
        S: Scalar;
      begin
        for K in reverse 0 .. D2 loop
          S := ScalZero;
          for I in IMax(0,K-P1.D) .. IMin(K,P2.D) loop
            S := S+P2.C(I)*P1.C(K-I);
          end loop;
          P2.C(K) := S;
        end loop;
        P2.D := D2;
      end;
      return;
    end if;                    -- done with numeric case

    declare                    -- now continue the easy way
      Pt: Taylor;
    begin
      Copy(P2,Pt);
      Prod(R,P1,Pt,P2,D);
    end;
  end Mult;

  procedure Prod(R: in Radius; P1,P2: in Taylor; P3: out Taylor; D: in Degree := Size) is
  begin
    if P1.D<0 or else P2.D<0 then
      P3.D := -1;
      return;
    end if;

    if Trunc then
      if D<0 then P3.D := -1; return; end if;
    elsif D<Dho then
      Prod(R,P1,P2,P3,Dho);
      DegConvert(R,D,P3);
      return;
    end if;

    declare
      D1: constant Natural := P1.D+P2.D;
      St,E: Scalar;
    begin
      P3.D := IMin(D1,D);
      for K in reverse 0 .. P3.D loop                  -- compute l.o. terms
        St := ScalZero;
        for I in IMax(0,K-P1.D) .. IMin(K,P2.D) loop
          St := St+P2.C(I)*P1.C(K-I);
        end loop;
        P3.C(K) := St;
      end loop;
      if Trunc then return; end if;                    -- done with numeric case

      if D1>D then                                     -- collect h.o. terms
        E := ScalZero;
        for K in reverse D+1 .. P2.D loop
          E := R*E+Ball(P2.C(K));
        end loop;
        St := (P1.C(0)+Ball(P1.E))*E;
        for K in reverse Imax(1,D+1-P1.D) .. Imin(D,P2.D) loop
          E := R*E+Ball(P2.C(K));
          St := St+P1.C(D+1-K)*E;
        end loop;
        if P1.D > D then
          E := R*E+Ball(P2.C(0))+Ball(P2.E);
          St := St+P1.C(D+1)*E;
          for K in D+2 .. P1.D loop
            E := R*E;
            St := St+P1.C(K)*E;
          end loop;
        end if;
        P3.C(D) := P3.C(D)+R*St;                       -- add them to P3.C(P3.D)
      end if;

      P3.E := P1.E*P2.E;                               -- lower order: err*err

      if P1.E /= ScalZero then                         -- lower order: err*coeff
        E := ScalZero;
        for K in reverse 0 .. Imin(Dho-1,P2.D) loop
          E := R*E+Abs(P2.C(K));
        end loop;
        P3.E := P3.E+P1.E*E;
        E := Ball(P1.E);
        for K in Dho .. Imin(P2.D,P3.D) loop
          P3.C(K) := P3.C(K)+E*P2.C(K);
        end loop;
      end if;

      if P2.E /= ScalZero then                         -- lower order: coeff*err
        E := ScalZero;
        for K in reverse 0 .. Imin(Dho-1,P1.D) loop
          E := R*E+Abs(P1.C(K));
        end loop;
        P3.E := P3.E+E*P2.E;
        E := Ball(P2.E);
        for K in Dho .. Imin(P1.D,P3.D) loop
          P3.C(K) := P3.C(K)+P1.C(K)*E;
        end loop;
      end if;
    end;
  end Prod;

  function Prod(R: Radius; P1,P2: Taylor; D: Degree := Size) return Taylor is
    Pt: Taylor;
  begin
    Prod(R,P1,P2,Pt,D);
    return Pt;
  end Prod;

  procedure MultAdd(R: in Radius; S: in Scalar; P1,P2: in Taylor; P3: in out Taylor; D: in Degree := Size; Monom1: in Boolean := False) is
  begin
    if P1.D<0 then                                     -- nothing to add
      DegConvert(R,D,P3);
      return;
    elsif Monom1 then                                  -- P1 is a monomial
      XMultAdd(R,P1.D,S*P1.C(P1.D),P2,P3,D);
      return;
    elsif P2.D<0 or else S=ScalZero then               -- nothing to add
      DegConvert(R,D,P3);
      return;
    end if;

    if Trunc then
      if D<0 then P3.D := -1; return; end if;
      if P1.D=0 then                                   -- P1 is a constant
        MultAdd(S*P1.C(0),P2,P3);
        DegConvert(R,D,P3);
        return;
      end if;
    elsif D<Dho then
      MultAdd(R,S,P1,P2,P3,Dho);
      DegConvert(R,D,P3);
      return;
    end if;

    declare
      D1: constant Natural := P1.D+P2.D;
      D2: constant Natural := IMin(D1,D);
      St,E: Scalar;
    begin
      if P3.D>D then
        DegConvert(R,D,P3);
      elsif P3.D<D2 then
        for K in P3.D+1 .. D2 loop
          P3.C(K) := ScalZero;
        end loop;
        if P3.D<0 then P3.E := ScalZero; end if;
        P3.D := D2;
      end if;

      for K in reverse 0 .. D2 loop                    -- compute l.o. terms
        St := ScalZero;
        for I in IMax(0,K-P1.D) .. IMin(K,P2.D) loop
          St := St+P2.C(I)*P1.C(K-I);
        end loop;
        P3.C(K) := P3.C(K)+S*St;
      end loop;
      if Trunc then return; end if;                    -- done with numeric case

      if D1>D then                                     -- collect h.o. terms
        E := ScalZero;
        for K in reverse D+1 .. P2.D loop
          E := R*E+Ball(P2.C(K));
        end loop;
        St := (P1.C(0)+Ball(P1.E))*E;
        for K in reverse Imax(1,D+1-P1.D) .. Imin(D,P2.D) loop
          E := R*E+Ball(P2.C(K));
          St := St+P1.C(D+1-K)*E;
        end loop;
        if P1.D > D then
          E := R*E+Ball(P2.C(0))+Ball(P2.E);
          St := St+P1.C(D+1)*E;
          for K in D+2 .. P1.D loop
            E := R*E;
            St := St+P1.C(K)*E;
          end loop;
        end if;
        P3.C(D) := P3.C(D)+R*S*St;                     -- add them to P3.C(P3.D)
      end if;

      P3.E := P3.E+Abs(S)*P1.E*P2.E;                   -- lower order: err*err

      if P1.E /= ScalZero then                         -- lower order: err*coeff
        E := ScalZero;
        for K in reverse 0 .. Imin(Dho-1,P2.D) loop
          E := R*E+Abs(P2.C(K));
        end loop;
        P3.E := P3.E+Abs(S)*P1.E*E;
        E := Ball(S*P1.E);
        for K in Dho .. Imin(P2.D,D2) loop
          P3.C(K) := P3.C(K)+E*P2.C(K);
        end loop;
      end if;

      if P2.E /= ScalZero then                         -- lower order: coeff*err
        E := ScalZero;
        for K in reverse 0 .. Imin(Dho-1,P1.D) loop
          E := R*E+Abs(P1.C(K));
        end loop;
        P3.E := P3.E+Abs(S)*E*P2.E;
        E := Ball(S*P2.E);
        for K in Dho .. Imin(P1.D,D2) loop
          P3.C(K) := P3.C(K)+P1.C(K)*E;
        end loop;
      end if;
    end;
  end MultAdd;

--- numeric stuff

  procedure Quot(P1,P2: in Taylor; P3: out Taylor; D: in Degree := Size) is
    R,S: Scalar;
  begin
    if not Trunc then
      Message("Taylors.Ops.Quot",Numeric_Only);
    end if;
    S := P2.C(0);
    if P2.D < 0 or else S = ScalZero then
      Message("Taylors.Ops.Quot",Division_By_Zero);
    end if;
    if D < 0 or else P1.D < 0 then
      P3.D := -1;
    else
      P3.D := D;
      P3.C(0) := P1.C(0)/S;
      for K in 1 .. D loop
        if K <= P1.D then
          R := P1.C(K);
        else
          R := ScalZero;
        end if;
        for I in 1 .. IMin(K,P2.D) loop
          R := R-P2.C(I)*P3.C(K-I);
        end loop;
        P3.C(K) := R/S;
      end loop;
    end if;
  end Quot;

  procedure Comp(P1,P2: in Taylor; P3: out Taylor) is
    R3: constant Radius := One;          --- arbitrary
  begin
    if not Trunc then
      P3.E := ScalZero;
      Message("Taylors.Ops.Comp",Numeric_Only);
    end if;
    if P1.D <= 0 then
      Copy(P1,P3);
    else
      Prod(P1.C(P1.D),P2,P3);
      P3.C(0) := P3.C(0)+P1.C(P1.D-1);
      for I in reverse 0 .. P1.D-2 loop
        Mult(R3,P2,P3);
        P3.C(0) := P3.C(0)+P1.C(I);
      end loop;
    end if;
  end Comp;

  function Comp(P1,P2: Taylor) return Taylor is
    Pt: Taylor;
  begin
    Comp(P1,P2,Pt);
    return Pt;
  end Comp;

  procedure InvComp(R,Err: in Radius; P1,P2: in Taylor; P3: in out Taylor) is
--- iterate p3 -> p3-[p2(p3)-p1]/p2'(p3)
    P,Pd,Pd2,Pd3: Taylor;
  begin
    if not Trunc then
      Message("Taylors.Ops.InvComp",Numeric_Only);
    end if;
    if P3.D < 0 then
      P3.D := 0;
      P3.C(0) := ScalZero;
      FindRoot(Err,P2,Val(R,ScalZero,P1),P3.C(0));
    end if;
    Der(R,R,P2,Pd2,True);
    for I in 1 .. 20 loop
      Comp(P2,P3,P);
      Sub(P1,P);
      Comp(Pd2,P3,Pd);
      Quot(P,Pd,Pd3);
      Sub(Pd3,P3);
      if Sup(Norm(R,Pd3)) <= Err then
        return;
      end if;
    end loop;
    Message("Taylors.Ops.InvComp",Giving_Up);
  end InvComp;

  function InvComp(R,Err: Radius; P1,P2: Taylor) return Taylor is
    Pt: Taylor;
  begin
    SetZero(Pt);
    InvComp(R,Err,P1,P2,Pt);
    return Pt;
  end InvComp;

  procedure FindRoot(Err: in Radius; P: in Taylor; S: in Scalar; R: in out Scalar) is
--- iterate r -> r+[s-p(r)]/p'(r)
    R1: constant Radius := One;       --- arbitrary
    Dr: Scalar;
  begin
    if not Trunc then
      Message("Taylors.Ops.FindRoot",Numeric_Only);
    end if;
    for I in 1 .. 20 loop
      Dr := (S-Val(R1,R,P))/DerVal(R1,R,P);
      R := R+Dr;
      if SupAbs(Dr) <= Err then
        return;
      end if;
    end loop;
    Message("Taylors.Ops.FindRoot",Giving_Up);
  end FindRoot;

end Taylors.Ops;