with Globals, Strings, Ints, Protected_Counters;
use Globals, Strings, Ints, Protected_Counters;

pragma Elaborate_All (Globals,Strings,Ints,Protected_Counters);

package body Polynoms2 is

  function IsSharp(P: Polynom2) return Boolean is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop
        if not IsSharp(P(I,J)) then return False; end if;
      end loop;
    end loop;
    return True;
  end IsSharp;

  function IsZero(P: Polynom2) return Boolean is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop
        if not IsZero(P(I,J)) then return False; end if;
      end loop;
    end loop;
    return True;
  end IsZero;

  function "="(P1,P2: Polynom2) return Boolean is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop
        if (P1(I,J) /= P2(I,J)) then return False; end if;
      end loop;
    end loop;
    return True;
  end "=";

  procedure SetZero(P: in out Triangular; Deg: in Positive) is
  begin
    for I in 0 .. Deg loop
      for J in 0 .. Deg-I loop SetZero(P(I,J)); end loop;
    end loop;
  end SetZero;

  procedure SetZero(P: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop SetZero(P(I,J)); end loop;
    end loop;
  end SetZero;

  procedure Copy(P1: in Polynom2; P2: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop Copy(P1(I,J),P2(I,J)); end loop;
    end loop;
  end Copy;

  procedure Swap(P1,P2: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop Swap(P1(I,J),P2(I,J)); end loop;
    end loop;
  end Swap;

  function Center0(P: Polynom2) return Boolean is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop
        if not Center0(P(I,J)) then return False; end if;
      end loop;
    end loop;
    return True;
  end Center0;

  function Contains0(P: Polynom2) return Logical is
    L: Logical := True;
    LN: Logical;
  begin
    for D in 1 .. PDeg loop
      for I in 0 .. D loop
        LN := Contains0(P(I,D-I));
        if LN=False then return False; end if;
        L := L and LN;
      end loop;
    end loop;
    return L;
  end Contains0;

  function Contains(P1,P2: Polynom2) return Logical is
    L: Logical := True;
    LN: Logical;
  begin
    for D in 1 .. PDeg loop
      for I in 0 .. D loop
        LN := Contains(P1(I,D-I),P2(I,D-I));
        if LN=False then return False; end if;
        L := L and LN;
      end loop;
    end loop;
    return L;
  end Contains;

  procedure ToErr(P: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop ToErr(P(I,J)); end loop;
    end loop;
  end ToErr;

  procedure ToErr(P1: in Polynom2; P2: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop ToErr(P1(I,J),P2(I,J)); end loop;
    end loop;
  end ToErr;

  procedure Center(P: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop Center(P(I,J)); end loop;
    end loop;
  end Center;

  procedure Center(P1: in Polynom2; P2: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop Center(P1(I,J),P2(I,J)); end loop;
    end loop;
  end Center;

  procedure ModCenter(P: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop ModCenter(P(I,J)); end loop;
    end loop;
  end ModCenter;

  procedure ModCenter(P1: in Polynom2; P2: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop ModCenter(P1(I,J),P2(I,J)); end loop;
    end loop;
  end ModCenter;

  procedure ErrMult(R: in Radius; P: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop ErrMult(R,P(I,J)); end loop;
    end loop;
  end ErrMult;

  procedure Union(P1: in Polynom2; P2: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop Union(P1(I,J),P2(I,J)); end loop;
    end loop;
  end Union;

  function QuasiDeg(P: Polynom2) return Natural is
  begin
    for D in reverse 1 .. PDeg loop
      for I in 0 .. D loop
        if not IsZero(P(I,D-I)) then return D; end if;
      end loop;
    end loop;
    return 0;
  end QuasiDeg;

  function QuasiDeg0(P: Polynom2) return Boolean is
  begin
    for D in reverse 1 .. PDeg loop
      for I in 0 .. D loop
        if not IsZero(P(I,D-I)) then return False; end if;
      end loop;
    end loop;
    return True;
  end QuasiDeg0;

  procedure Intersection(P1: in Polynom2; P2: in out Polynom2; Empty: out Logical) is
    E: Logical;
  begin
    Empty := False;
    for I in Power loop
      for J in 0 .. PDeg-I loop Intersection(P1(I,J),P2(I,J),E); end loop;
      Empty := Empty or E;
    end loop;
  end Intersection;

  procedure Neg(P: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop Neg(P(I,J)); end loop;
    end loop;
  end Neg;

  procedure Neg(P1: in Polynom2; P2: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop Neg(P1(I,J),P2(I,J)); end loop;
    end loop;
  end Neg;

  procedure Add(P1: in Triangular; P2: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop Add(P1(I,J),P2(I,J)); end loop;
    end loop;
  end Add;

  procedure Sum(P1,P2: in Polynom2; P3: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop
        Sum(P1(I,J),P2(I,J),P3(I,J));
      end loop;
    end loop;
  end Sum;

  procedure Sub(P1: in Polynom2; P2: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop Sub(P1(I,J),P2(I,J)); end loop;
    end loop;
  end Sub;

  procedure Diff(P1,P2: in Polynom2; P3: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop
        Diff(P1(I,J),P2(I,J),P3(I,J));
      end loop;
    end loop;
  end Diff;

  procedure Mult(R: in Flt; P: in out Polynom2; Tmp: in out Scalar) is
  begin
    if R=Zero then
      SetZero(P);
    elsif R=One then
      null;
    elsif R=NegOne then
      Neg(P);
    else
      for I in Power loop
        for J in 0 .. PDeg-I loop Mult(R,P(I,J),Tmp); end loop;
      end loop;
    end if;
  end Mult;

  procedure Prod(R: in Flt; P1: in Polynom2; P2: in out Polynom2) is
  begin
    if R=Zero then
      SetZero(P2);
    elsif R=One then
      Copy(P1,P2);
    elsif R=NegOne then
      Neg(P1,P2);
    else
      for I in Power loop
        for J in 0 .. PDeg-I loop
          Prod(R,P1(I,J),P2(I,J));
        end loop;
      end loop;
    end if;
  end Prod;

  procedure AddProd(R: in Flt; P1: in Polynom2; P2: in out Polynom2; Tmp: in out Scalar) is
  begin
    if R=Zero then
      null;
    elsif R=One then
      Add(P1,P2);
    elsif R=NegOne then
      Sub(P1,P2);
    else
      for I in Power loop
        for J in 0 .. PDeg-I loop
          AddProd(R,P1(I,J),P2(I,J),Tmp);
        end loop;
      end loop;
    end if;
  end AddProd;

  procedure Div(R: in Flt; P: in out Polynom2; Tmp: in out Scalar) is
  begin
    if R=One then
      null;
    elsif R=NegOne then
      Neg(P);
    else
      for I in Power loop
        for J in 0 .. PDeg-I loop Div(R,P(I,J),Tmp); end loop;
      end loop;
    end if;
  end Div;

  procedure Quot(P1: in Polynom2; R: in Flt; P2: in out Polynom2) is
  begin
    if R=One then
      Copy(P1,P2);
    elsif R=NegOne then
      Neg(P1,P2);
    else
      for I in Power loop
        for J in 0 .. PDeg-I loop Quot(P1(I,J),R,P2(I,J)); end loop;
      end loop;
    end if;
  end Quot;

  function IsReal(P: Polynom2) return Boolean is
  begin
    for D in Power loop
      for I in 0 .. D loop
        if not IsReal(P(I,D-I)) then return False; end if;
      end loop;
    end loop;
    return True;
  end IsReal;

  procedure Adjoint(P: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop
        Adjoint(P(I,J));
      end loop;
    end loop;
  end Adjoint;

  procedure Adjoint(P1: in Polynom2; P2: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop
        Adjoint(P1(I,J),P2(I,J));
      end loop;
    end loop;
  end Adjoint;

  procedure Real_Part(P: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop
        Real_Part(P(I,J));
      end loop;
    end loop;
  end Real_Part;

  procedure Imag_Part(P: in out Polynom2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop
        Imag_Part(P(I,J));
      end loop;
    end loop;
  end Imag_Part;

  procedure Mult(S: in Scalar; P: in out Polynom2) is
  begin
    if IsZero(S) then
      SetZero(P);
    else
      for I in Power loop
        for J in 0 .. PDeg-I loop Mult(S,P(I,J)); end loop;
      end loop;
    end if;
  end Mult;

  procedure Prod(S: in Scalar; P1: in Polynom2; P2: in out Polynom2) is
  begin
    if IsZero(S) then
      SetZero(P2);
    else
      for I in Power loop
        for J in 0 .. PDeg-I loop Prod(S,P1(I,J),P2(I,J)); end loop;
      end loop;
    end if;
  end Prod;

  procedure AddProd(S: in Scalar; P1: in Polynom2; P2: in out Polynom2) is
  begin
    if not IsZero(S) then
      declare
        Tmp: Scalar;
      begin
        for I in Power loop
          for J in 0 .. PDeg-I loop AddProd(S,P1(I,J),P2(I,J),Tmp); end loop;
        end loop;
      end;
    end if;
  end AddProd;

  procedure Add_Prod(P1,P2: in Polynom2; D1,D2,D3: in Natural; P3: in out Triangular) is
    Task_Error: Boolean := False;
    PC: Protected_Counter(0);

    procedure Prod_Loop is
      I3: Integer;
      Tmp: Scalar;
    begin
      Proper_Rounding;
      loop
        PC.Next(I3);
        exit when I3>D3;
        for J3 in 0 .. D3-I3 loop
          declare
            K: constant Integer := I3+J3-D2;
            S3: Scalar renames P3(I3,J3);
          begin
            for I1 in IMax(0,I3-D2) .. IMin(I3,D1) loop
              for J1 in IMax(0,K-I1) .. IMin(J3,D1-I1) loop
                AddProd(P1(I1,J1),P2(I3-I1,J3-J1),S3,Tmp);
              end loop;
            end loop;
          end;
        end loop;
      end loop;
    exception
      when others => Task_Error := True; raise;
    end Prod_Loop;

    task type Prod_Task_Type is end Prod_Task_Type;
    task body Prod_Task_Type is begin Prod_Loop; end Prod_Task_Type;

    NT: constant Natural := Reserve_Tasks(D3,PProd_Parallel);
  begin
    if NT=0 then
      Prod_Loop;
    else
      declare
        Prod_Task: array(1 .. NT) of Prod_Task_Type;
        pragma Warnings (Off,Prod_Task);
      begin
        Prod_Loop;
      end;
      Free_Tasks(NT);
      if Task_Error then raise Sorry with "Task_Error"; end if;
    end if;
  end Add_Prod;

  procedure Evaluate(P: in Polynom2; S1,S2: in Scalar; S3: in out Scalar; Deg: in Natural) is
  begin
    if Deg=0 then
      Copy(P(0,0),S3);
    else
      declare
        S4: Scalar;
      begin
        SetZero(S3);
        for I in reverse 0 .. Deg loop
          Copy(P(I,Deg-I),S4);
          for J in reverse 0 .. Deg-I-1 loop
            Mult(S2,S4);
            Add(P(I,J),S4);
          end loop;
          Mult(S1,S3);
          Add(S4,S3);
        end loop;
      end;
    end if;
  end Evaluate;

  procedure Evaluate(P: in Polynom2; S1,S2: in Scalar; S3: in out Scalar) is
  begin
    Evaluate(P,S1,S2,S3,QuasiDeg(P));
  end Evaluate;

  procedure Evaluate(P: in Polynom2; R1,R2: in Flt; S3: in out Scalar; Deg: in Natural) is
  begin
    if Deg=0 then
      Copy(P(0,0),S3);
    else
      declare
        S4: Scalar;
      begin
        SetZero(S3);
        for I in reverse 0 .. Deg loop
          Copy(P(I,Deg-I),S4);
          for J in reverse 0 .. Deg-I-1 loop
            Mult(R2,S4);
            Add(P(I,J),S4);
          end loop;
          Mult(R1,S3);
          Add(S4,S3);
        end loop;
      end;
    end if;
  end Evaluate;

  procedure Evaluate(P: in Polynom2; R1,R2: in Flt; S3: in out Scalar) is
  begin
    Evaluate(P,R1,R2,S3,QuasiDeg(P));
  end Evaluate;

  procedure Show1(N: in String; P: in Polynom2; Hide0: in Boolean := True) is
    W: constant Positive := 4;
    AllZero: Boolean := True;
    J: Integer;
  begin
    for D in 0 .. PDeg loop
      for I in 0 .. D loop
        J := D-I;
        if not IsZero(P(I,J)) then
          Show1(N & Strng(I,W) & Strng(J,W) & " ",P(I,J));
          AllZero := False;
        end if;
      end loop;
    end loop;
    if not Hide0 then
      if AllZero then Show0(N & "0"); end if;
      Show0;
    end if;
  end Show1;

  procedure Show2(N: in String; P1,P2: in Polynom2; Hide0: in Boolean := True) is
    W: constant Positive := 4;
    AllZero: Boolean := True;
    J: Integer;
  begin
    for D in 0 .. PDeg loop
      for I in 0 .. D loop
        J := D-I;
        if not (IsZero(P1(I,J)) and then IsZero(P2(I,J))) then
          Show2(N & Strng(I,W) & Strng(J,W) & " ",P1(I,J),P2(I,J));
          AllZero := False;
        end if;
      end loop;
    end loop;
    if not Hide0 then
      if AllZero then Show0(N & "0 0"); end if;
      Show0;
    end if;
  end Show2;

  procedure Num_Der1(P: in out Polynom2) is
  begin
    SetZero(P(0,PDeg));
    for I in 1 .. PDeg loop
      for J in 0 .. PDeg-I loop
        Prod(Flt(I),P(I,J),P(I-1,J));
      end loop;
      SetZero(P(I,PDeg-I));
    end loop;
  end Num_Der1;

  procedure Num_ValDer1(P: in Polynom2; S1,S2: in Scalar; D1: in out Scalar) is
    L: constant Integer := QuasiDeg(P);
  begin
    SetZero(D1);
    if L=0 then return; end if;
    declare
      S4,Tmp: Scalar;
    begin
      for I in reverse 1 .. L loop
        Prod(Flt(I),P(I,L-I),S4);
        for J in reverse 0 .. L-I-1 loop
          Mult(S2,S4);
          AddProd(Flt(I),P(I,J),S4,Tmp);
        end loop;
        Mult(S1,D1);
        Add(S4,D1);
      end loop;
    end;
  end Num_ValDer1;

  procedure Num_Der2(P: in out Polynom2) is
  begin
    SetZero(P(PDeg,0));
    for J in 1 .. PDeg loop
      for I in 0 .. PDeg-J loop
        Prod(Flt(J),P(I,J),P(I,J-1));
      end loop;
      SetZero(P(PDeg-J,J));
    end loop;
  end Num_Der2;

  procedure Num_ValDer2(P: in Polynom2; S1,S2: in Scalar; D2: in out Scalar) is
    L: constant Integer := QuasiDeg(P);
  begin
    SetZero(D2);
    if L=0 then return; end if;
    declare
      S4,Tmp: Scalar;
    begin
      for J in reverse 1 .. L loop
        Prod(Flt(J),P(L-J,J),S4);
        for I in reverse 0 .. L-J-1 loop
          Mult(S2,S4);
          AddProd(Flt(J),P(I,J),S4,Tmp);
        end loop;
        Mult(S1,D2);
        Add(S4,D2);
      end loop;
    end;
  end Num_ValDer2;

  procedure Num_MDer(P: in out Polynom2) is
    Tmp: Scalar renames P(0,0);
  begin
    for D in reverse 2 .. PDeg loop
      for I in 0 .. D loop
        Mult(Flt(D),P(I,D-I),Tmp);
      end loop;
    end loop;
    SetZero(P(0,0));
  end Num_MDer;

  procedure Inv_MDer(P: in out Polynom2; CheckConst: in Boolean := True) is
    Tmp: Scalar renames P(0,0);
  begin
    if CheckConst and then not IsZero(P(0,0)) then
      raise Undefined with "Matrices.Poly.Inv_MDer: nonzero constant";
    end if;
    for D in reverse 2 .. PDeg loop
      for I in 0 .. D loop
        Div(Flt(D),P(I,D-I),Tmp);
      end loop;
    end loop;
    SetZero(P(0,0));
  end Inv_MDer;

end Polynoms2;
