File : fouriers.ads

with Reps, Fouriers_Init;
use Reps, Fouriers_Init;
with Taylors;

generic

  Lmax, Wmax, Pdeg: in Positive;
  Dho: in Natural;
  type Scalar is private;

  with function IsNumeric(S: Scalar) return Boolean is <>;
  with function Scal(K: Integer) return Scalar is <>;
  with function Scal(R: Rep) return Scalar is <>;
  with function Ball(R: Rep) return Scalar is <>;
  with function Ball(S: Scalar) return Scalar is <>;
  with function Err0(S: Scalar) return Boolean is <>;
  with function "<"(S1,S2: Scalar) return Boolean is <>;
  with function IntFloor(S: Scalar) return Integer is <>;
  with function IntCeiling(S: Scalar) return Integer is <>;
  with function Approx(S: Scalar) return Rep is <>;
  with function Inf(S: Scalar) return Rep is <>;
  with function Sup(S: Scalar) return Rep is <>;
  with function SupAbs(S: Scalar) return Rep is <>;
  with function "abs"(S: Scalar) return Scalar is <>;
  with function Min(S1,S2: Scalar) return Scalar is <>;
  with function Max(S1,S2: Scalar) return Scalar is <>;
  with function Cap(R: Radius; S: Scalar) return Scalar is <>;
  with function "-"(S: Scalar) return Scalar is <>;

  with function Center(S: Scalar) return Scalar is <>;
  with function Widen(S: Scalar) return Scalar is <>;
  with procedure ResetCenter(S: in out Scalar) is <>;
  with function ResetCenter(S: Scalar) return Scalar is <>;
  with procedure Split(S: in Scalar; S0,SE: out Scalar) is <>;
  with procedure ErrMult(R: in Rep; S: in out Scalar) is <>;
  with function "+"(S,T: Scalar) return Scalar is <>;
  with function "-"(S,T: Scalar) return Scalar is <>;
  with function "*"(R: Rep; S: Scalar) return Scalar is <>;
  with function "/"(S: Scalar; R: Rep) return Scalar is <>;
  with function Sqr(S: Scalar) return Scalar is <>;
  with function "*"(S,T: Scalar) return Scalar is <>;
  with function Inv(S: Scalar) return Scalar is <>;
  with function "/"(S,T: Scalar) return Scalar is <>;
  with function "**"(S: Scalar; I: Integer) return Scalar is <>;

  with function Cosh(S: Scalar) return Scalar is <>;
  with function Sinh(S: Scalar) return Scalar is <>;
  with function Exp(S: Scalar) return Scalar is <>;
  with function Log(S: Scalar) return Scalar is <>;
  with function Lambda(M,N: Integer) return Scalar is <>;

package Fouriers is

  pragma Pure;

  Size1:  constant Positive := IntFloor(Lambda(Lmax,-Wmax)/Lambda(2,1));
  Size2:  constant Positive := IntFloor(Lambda(Wmax,-Lmax)/Lambda(2,1));
  SizeE:  constant Positive := NumErrs(Lmax,Wmax);

  subtype Freq1 is Integer range -Size1 .. Size1;
  subtype Freq2 is Integer range 0 .. Size2;
  subtype FreqE is Integer range 0 .. SizeE;
  subtype Indx3 is Integer range 0 .. Pdeg;

  package Components is new Taylors (Size => Pdeg, Dho => Dho, Scalar => Scalar);
  use Components;

  type Fourier is private;
  type FouVec is array (Integer range <>) of Fourier;

  function Deg(M,N: Integer; H: Fourier) return DegType;              -- Deg of H[M,N]
  function ErrDeg(I: Integer; H: Fourier) return DegType;             -- Deg of H_err[I]
  function Par(H: Fourier) return Parity;                             -- Parity of H  (0=even 1=odd)
  function Poly(H: Fourier;
                Check: in Boolean := False) return Boolean;           -- If true, coeffs are components
  procedure SetPoly(P: in Boolean; H: in out Fourier);                -- Declare H to be Poly or not
  function Rho(H: Fourier) return Weights;                            -- Weights defining domain of H
  procedure SetRho(R: in Weights; H: in out Fourier);                 -- Set Weights defining domain of H
  procedure SetZero(J: in Parity; R: in Weights; H: out Fourier);              -- H := 0

  procedure Component(M: in Freq1; N: in Freq2; H: in Fourier; P: out Taylor); -- P := H[M,N]
  function Component(M: Freq1; N: Freq2; H: Fourier) return Taylor;            -- H[M,N]
  function Component(M: Freq1; N: Freq2; K: Power; H: Fourier) return Scalar;  -- H[M,N,K]
  procedure ResetComponent(M: in Freq1; N: in Freq2; H: in out Fourier);       -- H[M,N] := 0
  procedure ErrComp(I: in FreqE; H: in Fourier; P: out Taylor);                -- P := H_err[I]
  function ErrComp(I: FreqE; H: Fourier) return Taylor;                        -- H_err[I]
  procedure ResetErrComp(I: in FreqE; H: in out Fourier);                      -- H_err[I] := 0
  procedure ResetErrComp(H: in out Fourier);                                   -- H_err := 0

  function NmaxDef return IntVec;                                     -- Standard cutoff for frequencies N
  function DmaxDef return Degrees;                                    -- Standard cutoff for degrees
  function DminusDef(Sig,Kap: Rep) return Degrees;                    -- Defines Iminus in RG_Ops

--- sobolev stuff

  type Sobolev is private;
  type SobVec is array (Integer range <>) of Sobolev;

  function Par(H: Sobolev) return Parity;                             -- Parity of H  (0=even 1=odd)
  function Poly(H: Sobolev) return Boolean;                           -- If true, H is Fourier polynomial
  function Rho(H: Sobolev) return Radius;                             -- Weights defining domain of H
  procedure SetRho(R: in Radius; H: in out Sobolev);                  -- Set Weights defining domain of H
  procedure SetZero(J: in Parity; R: in Radius; H: out Sobolev);      -- H := 0

  function Component(M: Freq1; N: Freq2; H: Sobolev) return Scalar;      -- H[M,N]
  procedure ResetComponent(M: in Freq1; N: in Freq2; H: in out Sobolev); -- H[M,N] := 0
  function ErrComp(I: FreqE; H: Sobolev) return Scalar;                  -- H_err[I]
  procedure ResetErrComp(I: in FreqE; H: in out Sobolev);                -- H_err[I] := 0
  procedure ResetErrComp(H: in out Sobolev);                             -- H_err := 0

private

  type Fourier is
    record
      P: Boolean;              --- if True, we have a Fourier-Taylor polynomial
      J: Parity;               --- even=0, odd=1
      R: Weights;              --- domain parameter
      C: TayMat(Freq1,Freq2);  --- coeffs are functions of z
      E: TayVec(FreqE);        --- these coeffs are also functions of q
    end record;

  type Sobolev is
    record
      J: Parity;               --- even=0, odd=1
      R: Radius;               --- domain parameter
      C: ScalMat(Freq1,Freq2); --- coefficients
      E: ScalVec(FreqE);       --- error terms
    end record;

  ScalZero: constant Scalar := Scal(Zero);
  ScalOne:  constant Scalar := Scal(One);
  Gold:     constant Scalar := Lambda(1,-1);
  InvGold:  constant Scalar := Lambda(0,-1);
  Trunc:   constant Boolean := IsNumeric(ScalZero);
  SizeE1: constant Positive := 2*Size1;
  SizeE2: constant Positive := 2*Size2+Size1;

  subtype FreqE1 is Integer range -SizeE1 .. SizeE1;
  subtype FreqE2 is Integer range 0 .. SizeE2;
  type IerrType is array (FreqE1,FreqE2) of FreqE;

  function EffSizeE(H: Fourier) return Integer;                            -- I of last nonzero H.E(I)
  function Ncut(M: Integer; L: Positive := Lmax;
                W: Positive := Wmax) return Integer;                       -- Standard cutoff for frequency N
  function Kcut(M,N: Integer; L: Positive := Lmax;
                W: Positive := Wmax; P: Positive := Pdeg) return DegType;  -- Standard cutoff for degree K
  function IerrDef return IerrType;                                        -- Indices I of H.E(I) for (M,N) errors

  function EffSizeE(H: Sobolev) return Integer;                            -- I of last nonzero H.E(I)

end Fouriers;