File : fouriers.adb

package body Fouriers is

  function Deg(M,N: Integer; H: Fourier) return DegType is
  begin
    return Deg(H.C(M,N));
  end Deg;

  function ErrDeg(I: Integer; H: Fourier) return DegType is
  begin
    if H.P then return -1; end if;
    return Deg(H.E(I));
  end ErrDeg;

  function Par(H: Fourier) return Parity is
  begin
    return H.J;
  end Par;

  function Poly(H: Fourier; Check: in Boolean := False) return Boolean is
  begin
    if Trunc then return True; end if;
    if not Check then return H.P; end if;
    if not H.P then return False; end if;
    for I in FreqE loop
      if Deg(H.E(I)) >= 0 then raise Constraint_Error; end if;
    end loop;
    return True;
  end Poly;

  procedure SetPoly(P: in Boolean; H: in out Fourier) is
  begin
    if Trunc then return; end if;
    H.P := P;
    if P then
      for M in Freq1 loop
        for N in Freq2 loop
          ResetErrComp(H.C(M,N));
        end loop;
      end loop;
      for I in FreqE loop
        SetZero(H.E(I));
      end loop;
    end if;
  end SetPoly;

  function Rho(H: Fourier) return Weights is
  begin
    return H.R;
  end Rho;

  procedure SetRho(R: in Weights; H: in out Fourier) is
  begin
    H.R := R;
  end SetRho;

  procedure SetZero(J: in Parity; R: in Weights; H: out Fourier) is
  begin
    H.P := True;
    H.J := J;
    H.R := R;
    for M in Freq1 loop
      for N in Freq2 loop
        SetZero(H.C(M,N));
      end loop;
    end loop;
    for I in FreqE loop
      SetZero(H.E(I));
    end loop;
  end SetZero;

  procedure Component(M: in Freq1; N: in Freq2; H: in Fourier; P: out Taylor) is
  begin
    Copy(H.C(M,N),P);
  end Component;

  function Component(M: Freq1; N: Freq2; H: Fourier) return Taylor is
  begin
    return H.C(M,N);
  end Component;

  function Component(M: Freq1; N: Freq2; K: Power; H: Fourier) return Scalar is
  begin
    return Component(K,H.C(M,N));
  end Component;

  procedure ResetComponent(M: in Freq1; N: in Freq2; H: in out Fourier) is
  begin
    SetZero(H.C(M,N));
  end ResetComponent;

  procedure ErrComp(I: in FreqE; H: in Fourier; P: out Taylor) is
  begin
    Copy(H.E(I),P);
  end ErrComp;

  function ErrComp(I: FreqE; H: Fourier) return Taylor is
  begin
    return H.E(I);
  end ErrComp;

  procedure ResetErrComp(I: in FreqE; H: in out Fourier) is
  begin
    SetZero(H.E(I));
  end ResetErrComp;

  procedure ResetErrComp(H: in out Fourier) is
  begin
    for I in FreqE loop
      SetZero(H.E(I));
    end loop;
  end ResetErrComp;

  function EffSizeE(H: Fourier) return Integer is
  begin
    for I in reverse FreqE loop
      if not IsZero(H.E(I)) then return I; end if;
    end loop;
    return -1;
  end EffSizeE;

  function Ncut(M: Integer; L: Positive := Lmax; W: Positive := Wmax) return Integer is
    N: Integer;
  begin
    N := Integer'Min(W-IntCeiling(Rep(M)*Inv(Gold)),IntFloor(Rep(L+M)*Gold));
    if N>Size2 then raise Constraint_Error; end if;
    return N;
  exception
    when Undefined => raise Constraint_Error; return N;
  end Ncut;

  function NmaxDef return IntVec is
    Nmax: IntVec(Freq1);
  begin
    for M in Freq1 loop
      Nmax(M) := Ncut(M);
    end loop;
    return Nmax;
  end NmaxDef;

  function Kcut(M,N: Integer; L: Positive := Lmax; W: Positive := Wmax; P: Positive := Pdeg) return DegType is
--- k+a*|omega*nu|<P+1, where a=(P+1-Dho)/W
--- k+b*|Omega*nu|<P+1, where b=(P+1-Dho)/L
    KW,KL,K: Integer;
    Nudotw,Nudotbigw: Scalar;
  begin
    K := -1;
    if M in Freq1 and then N in Nmin(M) .. Ncut(M,L,W) then
      Nudotw    := Abs(Rep( M)*InvGold+Scal(N));
      Nudotbigw := Abs(Scal(M)+Rep(-N)*InvGold);
      if IntFloor(Nudotw)<W and IntFloor(Nudotbigw)<L then
        KW := P-IntFloor(Rep(P+1-Dho)*Nudotw/Rep(W));
        KL := P-IntFloor(Rep(P+1-Dho)*Nudotbigw/Rep(L));
        K := Integer'Min(KW,KL);
        if K<Dho then raise Constraint_Error; end if;
      end if;
    end if;
    return K;
  exception
    when Undefined => raise Constraint_Error; return -1;
  end Kcut;

  function DmaxDef return Degrees is
    Dmax: Degrees(Freq1,Freq2);
  begin
    for M in Freq1 loop
      for N in Freq2 loop
        Dmax(M,N) := Kcut(M,N);
      end loop;
    end loop;
    return Dmax;
  end DmaxDef;

  function IerrDef return IerrType is
    Nudotw,Nudotbigw: Scalar;
    Ierr: IerrType;
  begin
    for M in FreqE1 loop
      for N in FreqE2 loop
        Nudotw    := Abs(Rep( M)*InvGold+Scal(N));
        Nudotbigw := Abs(Scal(M)+Rep(-N)*InvGold);
        Ierr(M,N) := Integer'Min(SizeE,Integer'Max(IntFloor(Nudotw),IntFloor(Nudotbigw)));
      end loop;
    end loop;
    return Ierr;
  exception
    when Undefined => raise Constraint_Error; return Ierr;
  end IerrDef;

  function DminusDef(Sig,Kap: Rep) return Degrees is
    Nudotw,Nudotbigw: Scalar;
    Dminus: Degrees(Freq1,Freq2);
  begin
    for M in Freq1 loop
      for N in Freq2 loop
        Nudotw    := Abs(Rep( M)*InvGold+Scal(N));
        Nudotbigw := Abs(Scal(M)+Rep(-N)*InvGold);
        if Sig*Nudotbigw < Nudotw then
          Dminus(M,N) := IntFloor(Nudotw/Kap);
        else
          Dminus(M,N) := -1;
        end if;
      end loop;
    end loop;
    return Dminus;
  exception
    when Undefined => raise Constraint_Error; return Dminus;
  end DminusDef;

end Fouriers;