with Globals;
use Globals;

pragma Elaborate_All (Globals);

package body Vectors.More is

  procedure ErrMult(R: in Radius; V: in out Vector) is
  begin
    for N in V'Range loop ErrMult(R,V(N)); end loop;
  end ErrMult;

  procedure Intersection(V1: in Vector; V2: in out Vector; Empty: out Logical) is
    E: Logical;
  begin
    if CheckRange and then ((V1'First /= V2'First) or (V1'Last /= V2'Last)) then
      raise Dimension_Error with "Vectors.More.Intersection error";
    end if;
    Empty := False;
    for N in V1'Range loop
      Intersection(V1(N),V2(N),E);
      Empty := Empty or E;
    end loop;
  end Intersection;

  procedure Min(V: in Vector; S: in out Component) is
  begin
    for N in V'Range loop Min(V(N),S); end loop;
  end Min;

  procedure Max(V: in Vector; S: in out Component) is
  begin
    for N in V'Range loop Max(V(N),S); end loop;
  end Max;

  function Sup(V: Vector) return Flt is
    R: Flt := Sup(V(V'First));
  begin
    for N in V'First+1 .. V'Last loop
      R := RMax(Sup(V(N)),R);
    end loop;
    return R;
  end Sup;

  function Inf(V: Vector) return Flt is
    R: Flt := Inf(V(V'First));
  begin
    for N in V'First+1 .. V'Last loop
      R := RMin(Inf(V(N)),R);
    end loop;
    return R;
  end Inf;

  procedure Mult(R: in Flt; V: in out Vector; Tmp: in out Component) is
  begin
    if R=Zero then
      for N in V'Range loop SetZero(V(N)); end loop;
    elsif R=One then
      null;
    elsif R=NegOne then
      for N in V'Range loop Neg(V(N)); end loop;
    else
      for N in V'Range loop Mult(R,V(N),Tmp); end loop;
    end if;
  end Mult;

  procedure Prod(R: in Flt; V1: in Vector; V2: in out Vector) is
  begin
    if CheckRange and then ((V1'First /= V2'First) or (V1'Last /= V2'Last)) then
      raise Dimension_Error with "Vectors.More.Prod error";
    end if;
    if R=Zero then
      for N in V1'Range loop SetZero(V2(N)); end loop;
    elsif R=One then
      for N in V1'Range loop Copy(V1(N),V2(N)); end loop;
    elsif R=NegOne then
      for N in V1'Range loop Neg(V1(N),V2(N)); end loop;
    else
      for N in V1'Range loop Prod(R,V1(N),V2(N)); end loop;
    end if;
  end Prod;

  procedure AddProd(R: in Flt; V1: in Vector; V2: in out Vector; Tmp: in out Component) is
  begin
    if CheckRange and then ((V1'First /= V2'First) or (V1'Last /= V2'Last)) then
      raise Dimension_Error with "Vectors.More.AddProd error";
    end if;
    if R=Zero then
      null;
    elsif R=One then
      for N in V1'Range loop Add(V1(N),V2(N)); end loop;
    elsif R=NegOne then
      for N in V1'Range loop Sub(V1(N),V2(N)); end loop;
    else
      for N in V1'Range loop AddProd(R,V1(N),V2(N),Tmp); end loop;
    end if;
  end AddProd;

  procedure Div(R: in Flt; V: in out Vector; Tmp: in out Component) is
  begin
    if R=One then
      null;
    elsif R=NegOne then
      for N in V'Range loop Neg(V(N)); end loop;
    else
      for N in V'Range loop Div(R,V(N),Tmp); end loop;
    end if;
  end Div;

  procedure Quot(V1: in Vector; R: in Flt; V2: in out Vector) is
  begin
    if CheckRange and then ((V1'First /= V2'First) or (V1'Last /= V2'Last)) then
      raise Dimension_Error with "Vectors.More.Quot error";
    end if;
    if R=One then
      for N in V1'Range loop Copy(V1(N),V2(N)); end loop;
    elsif R=NegOne then
      for N in V1'Range loop Neg(V1(N),V2(N)); end loop;
    else
      for N in V1'Range loop Quot(V1(N),R,V2(N)); end loop;
    end if;
  end Quot;

  procedure Evaluate(P: in Polynom1; R: in Flt; C: in out Component; Deg: in Natural) is
  begin
    if Deg=0 or else R=Zero then
      Copy(P(0),C);
    else
      Copy(P(Deg),C);
      for I in reverse 0 .. Deg-1 loop
        Mult(R,C);
        Add(P(I),C);
      end loop;
    end if;
  end Evaluate;

  procedure Evaluate(P: in Polynom1; R: in Flt; C: in out Component) is
  begin
    Evaluate(P,R,C,EffLast(P));
  end Evaluate;

  function Evaluate(P: Polynom1; R: Flt; Deg: Natural) return Component is
    C: Component;
  begin
    Evaluate(P,R,C,Deg);
    return C;
  end Evaluate;

  function Evaluate(P: Polynom1; R: Flt) return Component is
    C: Component;
  begin
    Evaluate(P,R,C,EffLast(P));
    return C;
  end Evaluate;

--- derivatives etc

  procedure Num_Der(P: in out Polynom1) is
  begin
    for N in 1 .. P'Last loop
      Prod(Flt(N),P(N),P(N-1));
    end loop;
    SetZero(P(P'Last));
  end Num_Der;

  procedure Num_MDer(P: in out Polynom1) is
    Tmp: Component renames P(0);
  begin
    for N in reverse 2 .. P'Last loop
      Mult(Flt(N),P(N),Tmp);
    end loop;
    SetZero(P(0));
  end Num_MDer;

  procedure Inv_MDer(P: in out Polynom1; CheckConst: in Boolean := True) is
    Tmp: Component renames P(0);
  begin
    if CheckConst and then not IsZero(P(0)) then
      raise Undefined with "Vectors.More.Inv_MDer: nonzero constant";
    end if;
    for N in reverse 2 .. P'Last loop
      Div(Flt(N),P(N),Tmp);
    end loop;
    SetZero(P(0));
  end Inv_MDer;

end Vectors.More;
