with Ada.Text_IO,MultiScalars, Vectors;
use Ada.Text_IO;

pragma Elaborate_All (Ada.Text_IO,MultiScalars,Vectors);

generic

  type Scalar is private;

  --- for Vectors
  with function "="(C1,C2: Scalar) return Boolean is <>;
  with function IsZero(C: Scalar) return Boolean is  <>;
  with function IsReal(C: Scalar) return Boolean is  <>;
  with procedure SetZero(C: in out Scalar) is <>;
  with procedure Copy(C1: in Scalar; C2: in out Scalar) is <>;
  with procedure Swap(C1,C2: in out Scalar) is <>;
  with procedure Neg(C: in out Scalar) is <>;
  with procedure Neg(C1: in Scalar; C2: in out Scalar) is <>;
  with procedure Add(C1: in Scalar; C2: in out Scalar) is <>;
  with procedure Sum(C1,C2: in Scalar; C3: in out Scalar) is <>;
  with procedure Sub(C1: in Scalar; C2: in out Scalar) is <>;
  with procedure Diff(C1,C2: in Scalar; C3: in out Scalar) is <>;
  with procedure Real_Part(C: in out Scalar) is <>;
  with procedure Imag_Part(C: in out Scalar) is <>;
  with procedure Adjoint(C: in out Scalar) is <>;
  with procedure Adjoint(C1: in Scalar; C2: in out Scalar) is <>;
  with procedure Show1(N: in String; C: in Scalar; Whatever: in Boolean := True) is <>;
  with procedure Show2(N: in String; C1,C2: in Scalar; Whatever: in Boolean := True) is <>;
  with procedure Put(F: in File_Type; C: in Scalar; Decimal: in Boolean := False) is <>;
  with procedure Get(F: in File_Type; C: in out Scalar; Decimal: in Boolean := False) is <>;

  --- more
  with procedure Assign(I: in Integer; S: in out Scalar) is <>;
  with procedure Mult(S1: in Scalar; S2: in out Scalar) is <>;
  with procedure Prod(S1,S2: in Scalar; S3: in out Scalar) is <>;
  with procedure AddProd(S1,S2: in Scalar; S3: in out Scalar) is <>;

  with package SV is new Vectors (Component => Scalar);

package Polynoms is

  subtype Vector is SV.Vector;
  type Polynom1 is new Vector; --- Natural range

  function EffDeg(P: Polynom1) return Natural;
  function EffDeg0(P: Polynom1) return Boolean;
  function IsEven(C: Polynom1) return Boolean;
  function IsOdd(C: Polynom1) return Boolean;
  procedure Copy(P1: in Polynom1; P2: in out Polynom1);
  procedure Mult(S: in Scalar; P: in out Polynom1);
  procedure Prod(S: in Scalar; P1: in Polynom1; P2: in out Polynom1);
  procedure AddProd(S: in Scalar; P1: in Polynom1; P2: in out Polynom1);
  procedure Evaluate(P: in Polynom1; S: in Scalar; C: in out Scalar; Deg: in Natural); -- value of P at S
  procedure Evaluate(P: in Polynom1; S: in Scalar; C: in out Scalar);                  -- value of P at S
  function Evaluate(P: Polynom1; S: Scalar; Deg: Natural) return Scalar;               -- value of P at S
  function Evaluate(P: Polynom1; S: Scalar) return Scalar;                             -- value of P at S
  procedure Evaluate(P: in Polynom1; Sr,Si: in Scalar; Cr,Ci: in out Scalar);          -- value of P at Sr+i*Si
  procedure EvalDer(P: in Polynom1; S: in Scalar; C: in out Scalar; Deg: in Natural);  -- value of P' at S
  procedure EvalDer(P: in Polynom1; S: in Scalar; C: in out Scalar);                   -- value of P' at S

  procedure Scale(S: in Scalar; P: in out Polynom1);                           -- P(x) := P(S*x)
  procedure Scale(P1: in Polynom1; S: in Scalar; P2: in out Polynom1);         -- P2(x) := P1(S*x)
  procedure Translate(P1: in Polynom1; S: in Scalar; P2,Tmp: in out Polynom1); -- P2(x) := P1(x+S)
  procedure Mult1(S: in Scalar; P: in out Polynom1);                           -- P(x) := (S-x)*P(x)
  procedure Prod1(S: in Scalar; P1: in Polynom1; P2: in out Polynom1);         -- P2(x) := (S-x)*P1(x)
  procedure Prod(P1,P2: in Polynom1; P3: in out Polynom1);                     -- P3 := P1*P2

  procedure LMultiply(W: in Scalar; P: in out Polynom1);                      -- P(x) := (x-W)*P(x)
  procedure LMultiply(A,B: in Scalar; P: in out Polynom1);                    -- P(x) := (A*x+B)*P(x)
  procedure QMultiply(U,V,C,D: in Scalar; P: in out Polynom1);                -- P(x) := (x*x+U*x+V)*P(x)+(C*x+D)
  procedure QProduct(U,V,C,D: in Scalar; Q: in Polynom1; P: in out Polynom1); -- P(x) := (x*x+U*x+V)*Q(x)+(C*x+D)
  procedure QMultiply(U,V: in Vector; A,B: in Scalar; P: in out Polynom1);    -- P(x) := (A*x+B)*Prod_k(x*x+U(K)*x+V(K))
  procedure LDivision(P: in out Polynom1; W: in Scalar; R: in out Scalar);                     -- divivion of P(x) by by (x-W)
  procedure LDivision(P: in Polynom1; W: in Scalar; Q: in out Polynom1; R: in out Scalar);     -- divivion of P(x) by by (x-W)
  procedure QDivision(P: in Polynom1; U,V: in Scalar; Q: in out Polynom1; C,D: in out Scalar); -- division of P(x) by (x*x+U*x+V)
  procedure QDivision(P: in Polynom1; U,V: in Scalar; Q: in out Polynom1);                     -- division of P(x) by (x*x+U*x+V)
  procedure Roots2Poly(W: in Vector; P: in out Polynom1);                     -- P with real roots W(I)
  procedure Roots2Poly(Wr,Wi: in Vector; P: in out Polynom1);                 -- P with roots Wr(I) plusminus i*Wi(I)

  subtype Cheby is Vector;
  procedure Chebyshev(N: in Natural; P: in out Polynom1);      -- N-th Chebyshev polynomial
  procedure Clenshaw(C: in Cheby; P: in out Polynom1);         -- convert finite Cheby series to polynomial
  procedure SymPolyCheby(P: in Polynom1; C: in out Cheby);     -- Cheby assiciated with P
  procedure SymCosPoly(P: in Polynom1; Q: in out Polynom1);    -- roots of Q are (r+1/r)/2

  type Polynom2 is new Vector; --- index = J+D*(D+1)/2   where D=I+J

  function MaxDeg(P: Polynom2) return Natural;
  function EffDeg(P: Polynom2) return Natural;

  procedure SetCoeff(I,J: in Natural; C: in Scalar; P: in out Polynom2);
  procedure Coefficient(P: in Polynom2; I,J: in Natural; C: in out Scalar);
  procedure Show1(N: in String; P: in Polynom2; Hide0: in Boolean := True);
  procedure Show2(N: in String; P1,P2: in Polynom2; Hide0: in Boolean := True);
  procedure Evaluate(P: in Polynom2; S1,S2: in Scalar; C: in out Scalar);              -- value of P at (S1,S2)

  package MuSc is new MultiScalars (Scalar => Scalar);

private

  use MuSc;

  procedure Bairstow0(P: in Polynom1; Deg: in Integer; S: in Vector2; Q: in out Polynom1; FS: in out Vector2);  -- QDivide map F
  procedure Bairstow1(Q: in Polynom1; Deg: in Integer; S: in Vector2; R: in out Polynom1; DFS: in out Matrix2); -- derivative DF

end Polynoms;
