File : taylors-ops.ads

with Reps;
use Reps;

generic

package Taylors.Ops is

  procedure ComponentSplit(K: in Power; P: in Taylor;
                           P0: out Taylor; E: out Scalar);             -- Split P[K] into P0[0]+E
  procedure AddComponent(D: in Power; S: in Scalar; P: in out Taylor); -- P[D] := P[D] + S
  procedure MultComponent(D: in Power; S: in Scalar; P: in out Taylor);-- P[D] := P[D] + S
  procedure AddErrComp(S: in Scalar; P: in out Taylor);                -- add |S| to general error
  procedure MultErrComp(S: in Scalar; P: in out Taylor);               -- add |S| to general error
  function Coeff(R: Radius; K: Power; P: Taylor) return Scalar;        -- P_K
  procedure ResetCenter(P: in out Taylor;
                   Kmin: in Natural := 0; Kmax: in Integer := Size);   -- Reset centers of P.C(Kmin .. Kmax)
  function Val(R: Radius; Z: Scalar; P: Taylor) return Scalar;         -- P(Z)
  function Norm(R: Radius; P: Taylor) return Scalar;                   -- Norm ||P||_R
  procedure ErrNorm(R: in Radius; Pw,P1: in Taylor;
                    P2: out Taylor; E: out Scalar);                    -- Used by Norm below, and ...

  procedure Neg(P1: in Taylor; P2: out Taylor);                        -- P2 := -P1
  function "-"(P: Taylor) return Taylor;                               -- -P
  procedure Mult(S: in Scalar; P: in out Taylor);                      -- P := S*P
  procedure Prod(S: in Scalar; P1: in Taylor; P2: out Taylor);         -- P2 := S*P1
  function "*"(S: Scalar; P: Taylor) return Taylor;                    -- S*P
  procedure Add(P1: in Taylor; P2: in out Taylor);                     -- P2 := P2+P1
  procedure Sum(P1,P2: in Taylor; P3: out Taylor);                     -- P3 := P1+P2
  function "+"(P1,P2: Taylor) return Taylor;                           -- P1+P2
  procedure Sub(P1: in Taylor; P2: in out Taylor);                     -- P2 := P2-P1
  procedure Diff(P1,P2: in Taylor; P3: out Taylor);                    -- P3 := P1-P2
  function "-"(P1,P2: Taylor) return Taylor;                           -- P1-P2
  procedure MultAdd(S: in Scalar; P1: in Taylor; P2: in out Taylor);   -- P2 := P2+S*P1
  procedure Scale(S: in Scalar; P1: in Taylor; P2: out Taylor);        -- P2 := P1(S.)/S
  procedure Scale(S: in Scalar; P: in out Taylor);                     -- P := P(S.)/S
  function Scale(S: Scalar; P: Taylor) return Taylor;                  -- Scale := P(S.)/S
  procedure ScaleRem1(Eps: in Scalar; P1: in Taylor; P2: out Taylor;
                      IsPoly: in Boolean := False);                    -- P2 := (1+Eps)*P1(./(1+Eps))-P1
  procedure ScaleRem2(Eps: in Scalar; P1: in Taylor; P2: out Taylor;   -- P2 := (1+Eps)*P1(./(1+Eps))-P1
                      IsPoly: in Boolean := False);                    --       -Eps*P1+Eps*Z*P1'

  procedure Fmax1(R: in Scalar; M: out Natural; F: out Scalar);        -- used by Der
  procedure Der(R1,R2: in Radius; P1: in Taylor; P2: out Taylor;
                IsPoly: in Boolean := False);                          -- P2 := P1'
  function Der(R1,R2: Radius; P: Taylor;
               IsPoly: in Boolean := False) return Taylor;             -- Der := P'
  procedure DerMult(R,Fmax: in Radius; S: in Scalar;
                    P: in out Taylor; Plain: in Boolean := False);     -- S*P' and projection
  procedure DerProd(R,Fmax: in Radius; S: in Scalar; P1: in Taylor;
                    P2: out Taylor; Plain: in Boolean := False);       -- P2 := S*P1' and projection
  procedure ZxDer(P1: in Taylor; P2: out Taylor;
                  IsPoly: in Boolean := False);                        -- P2 := (.)*P1'
  procedure DerDer(P: in out Taylor; IsPoly: in Boolean := False);     -- P := P"

  procedure DegConvert(R: in Radius; D: in Degree; P: in out Taylor);  -- Enlarge P by reducing Deg
  procedure Mult(R: in Radius; P1: in Taylor; P2: in out Taylor;
                 D: in Degree := Size);                                -- P2 := P2*P1
  procedure Prod(R: in Radius; P1,P2: in Taylor; P3: out Taylor;
                 D: in Degree := Size);                                -- P3 := P1*P2
  function Prod(R: Radius; P1,P2: Taylor;
                D: Degree := Size) return Taylor;                      -- P1*P2
  procedure MultAdd(R: in Radius; S: in Scalar; P1,P2: in Taylor;
                    P3: in out Taylor; D: in Degree := Size;
                    Monom1: in Boolean := False);                      -- P3 := P3+S*P1*P2

--- the following are meant to be used only with Numeric=True

  procedure Quot(P1,P2: in Taylor; P3: out Taylor; D: in Degree := Size);  -- P3 := P1/P2
  procedure Comp(P1,P2: in Taylor; P3: out Taylor);                    -- P3 := P1 o P2
  function Comp(P1,P2: Taylor) return Taylor;                          -- P1 o P2
  procedure InvComp(R,Err: in Radius; P1,P2: in Taylor;
                    P3: in out Taylor);                                -- P3 solves P2(P3)=P1
  function InvComp(R,Err: Radius; P1,P2: Taylor) return Taylor;        -- InvComp solves P2(InvComp)=P1
  function DerVal(R: Radius; Z: Scalar; P: Taylor) return Scalar;      -- P'(Z)
  procedure FindRoot(Err: in Radius; P: in Taylor; S: in Scalar;
                     R: in out Scalar);                                -- R solves P(R)=S

  pragma Inline (AddComponent,MultComponent,AddErrComp,MultErrComp,Coeff);

private

  procedure XMultAdd(R: in Radius; K: in Power; S: in Scalar;
             P1: in Taylor; P2: in out Taylor; D: in Degree := Size); -- P2(Z) := P2(Z)+S*(Z^K)*P1(Z)

end Taylors.Ops;