File : taylors.adb

package body Taylors is

  function Deg(P: Taylor) return Degree is
  begin
    return P.D;
  end Deg;

  function EffDeg(P: Taylor) return Degree is
  begin
    if P.D<0 then return -1; end if;
    for K in reverse 0 .. P.D loop
      if P.C(K) /= ScalZero then
        return K;
      end if;
    end loop;
    if Trunc or else P.E=ScalZero then return -1; end if;
    return 0;
  end EffDeg;

  procedure SetDeg(D: in Degree; P: in out Taylor) is
  begin
    if D<P.D then P.D := D; end if;
  end SetDeg;

  procedure SetZero(P: out Taylor) is
  begin
    P.D := -1;
  end SetZero;

  function IsZero(P: Taylor) return Boolean is
  begin
    if P.D<0 then return True; end if;
    for K in reverse 0 .. P.D loop
      if P.C(K) /= ScalZero then
        return False;
      end if;
    end loop;
    return Trunc or else P.E=ScalZero;
  end IsZero;

  procedure Cleanup(P: in out Taylor) is
  begin
    if P.D<0 then return; end if;
    for K in reverse 0 .. P.D loop
      exit when P.C(K) /= ScalZero;
      P.D := K-1;
    end loop;
    if Trunc or else P.D >= 0 or else P.E=ScalZero then return; end if;
    P.D := 0;
  end Cleanup;

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

  function Component(K: Power; P: Taylor) return Scalar is
  begin
    if K > P.D then
      return ScalZero;
    end if;
    return P.C(K);
  end Component;

  procedure Component(K: in Power; P1: in Taylor; P2: out Taylor) is
  begin
    if K>P1.D then
      P2.D := -1;
    else
      P2.D := 0;
      P2.C(0) := P1.C(K);
      P2.E := ScalZero;
    end if;
  end Component;

  procedure SetComponent(D: in Power; S: in Scalar; P: in out Taylor) is
  begin
    if D>P.D then
      if S=ScalZero then return; end if;
      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;
      return;
    end if;
    if D /= P.D or else S /= ScalZero then
      P.C(D) := S;
      return;
    end if;
    if Trunc or else D>0 or else P.E=ScalZero then
      P.D := D-1;
      return;
    end if;
    P.C(0) := S;                  -- when only P.E remains nonzero
  end SetComponent;

  procedure ResetComponent(P: in out Taylor; Kmin: in Natural := 0; Kmax: in Integer := Size) is
  begin
    if Kmax<P.D then
      for K in Kmin .. Kmax loop
        P.C(K) := ScalZero;
      end loop;
    elsif Kmin<=P.D then
      if Trunc or else Kmin>0 or else P.E=ScalZero then
        P.D := Kmin-1;
      else
        P.C(0) := ScalZero;    -- when only P.E remains nonzero
      end if;
    end if;
  end ResetComponent;

  function ErrComp(P: Taylor) return Scalar is
  begin
    if Trunc or else P.D<0 then
      return ScalZero;
    end if;
    return P.E;
  end ErrComp;

  procedure SetErrComp(S: in Scalar; P: in out Taylor) is
  begin
    if Trunc then return; end if;
    if P.D<0 then
      if S=ScalZero then return; end if;
      P.D := 0;
      P.C(0) := ScalZero;        -- when only P.E remains
    end if;
    P.E := Abs(S);
  end SetErrComp;

  procedure ResetErrComp(P: in out Taylor) is
  begin
    if Trunc then return; end if;
    if P.D=0 and then P.C(0)=ScalZero then P.D := -1; return; end if;
    P.E := ScalZero;
  end ResetErrComp;

end Taylors;