
package body Generic_DFT is

  function Even(P: FPoly) return Boolean is
  begin
     return True;
  end Even;

  procedure Proj_Even(P: in out FPoly) is
  begin
    for I in reverse 1-L .. -1 loop
      SetZero(P(I));
    end loop;
  end Proj_Even;

  procedure Val(P: in FPoly; X: in Integer; V,Tmp: in out Scalar) is
    NX: Integer;
  begin
    Copy(P(L),V);
    if (X mod 2)=1 then Neg(V); end if;
    for N in reverse NRange loop
      NX := (N*X) mod T;
      AddProd(P( N),ValCos(NX),V,Tmp);
    end loop;
    Add(P(0),V);
  end Val;

  procedure Val(P: in FPoly; V: in out Values) is
    Tmp: Scalar;
  begin
    for X in XRange loop Val(P,X,V(X),Tmp); end loop;
  end Val;

  procedure CosCoeff(N: in Natural; V: in Values; S,Tmp: in out Scalar) is
  begin
    Assign(Zero,S);
    for X in XRange loop
      AddProd(V(X),ValCos((N*X) mod T),S,Tmp);
    end loop;
    if N=0 or else N=L then
      Div(Flt(2*L),S,Tmp);
    else
      Div(Flt(L),S,Tmp);
    end if;
  end CosCoeff;

  procedure SinCoeff(N: in Natural; V: in Values; S,Tmp: in out Scalar) is
  begin
    Assign(Zero,S);
    for X in XRange loop
      AddProd(V(X),ValSin((N*X) mod T),S,Tmp);
    end loop;
    Div(Flt(L),S,Tmp);
  end SinCoeff;

  procedure DCT(V: in Values; P: in out FPoly) is
    Tmp: Scalar;
  begin
    CosCoeff(0,V,P(0),Tmp);
    for N in NRange loop
      CosCoeff(N,V,P(N),Tmp);
    end loop;
    CosCoeff(L,V,P(L),Tmp);
  end DCT;

  procedure DFT(V: in Values; P: in out FPoly) is
    Tmp0,Tmp1: Scalar;
    task Run_SinCoeff is end Run_SinCoeff;
    task body Run_SinCoeff is
    begin
      for N in NRange loop SinCoeff(N,V,P(-N),Tmp1); end loop;
    end Run_SinCoeff;
  begin
    for N in 0 .. L loop CosCoeff(N,V,P(N),Tmp0); end loop;
  end DFT;

  procedure F1(P: in out FPoly) is
    V1,V2: Values;
  begin
    Val(P,V1);
    for X in XRange loop
      F(V1(X),V2(X));
    end loop;
    if Even(P) then
      DCT(V2,P);
    else
      DFT(V2,P);
    end if;
  end F1;

  procedure F2(P1: in FPoly; P2: in out FPoly) is
    V1,V2: Values;
  begin
    Val(P1,V1);
    for X in XRange loop
      F(V1(X),V2(X));
    end loop;
    if Even(P1) then
      DCT(V2,P2);
    else
      DFT(V2,P2);
    end if;
  end F2;

begin
  SetZero(STmp);
  ArcCos(STmp,PiOverL);       -- Pi/2
  Mult(Flt(2),PiOverL,CTmp);  -- Pi
  Div(Flt(L),PiOverL,CTmp);   -- Pi/L

  for X in 1 .. T-1 loop
    Copy(PiOverL,CTmp);
    Mult(Flt(X),CTmp,STmp);
    Cos(CTmp,ValCos(X));     -- Cos(X*Pi/L)
    Sin(CTmp,ValSin(X));     -- Sin(X*Pi/L)
  end loop;
  Assign(One,ValCos(0));
  Assign(Zero,ValSin(0));

end Generic_DFT;
