File : fouriers-ops.ads


with Taylors.Ops;

pragma Elaborate_All (Taylors.Ops);

generic

package Fouriers.Ops is

  Nmax: IntVec(Freq1) := NmaxDef;
  Dmax: Degrees(Freq1,Freq2) := DmaxDef;

  type FouVecPtr is access FouVec;

  package Components_Ops is new Components.Ops;

  function DmaxCount return Integer;                                  -- Number of coefficients
  function ErrCount return Integer;                                   -- Number of error terms
  function InRange(M,N: Integer; K: Natural := 0;
                   D: Degrees := Dmax) return Boolean;                -- InRange := (M and N are in range)

  procedure Component(K: Natural; H1: in Fourier; H2: out Fourier);   -- H2[M,N,0] := H1[M,N,K]  (roughly)
  procedure SetComponent(M: in Freq1; N: in Freq2; K: in Natural;
                     S: in Scalar; H: in out Fourier);                -- H[M,N,K] := S
  procedure AddComponent(M: in Freq1; N: in Freq2; K: in Natural;
                     S: in Scalar; H: in out Fourier);                -- Add S to H[M,N,K]
  procedure SetErrComp(I: in FreqE; K: in Natural; S: in Scalar;
                       H: in out Fourier);                            -- H_err[I,K] := Ball(S)
  procedure AddErrComp(I: in FreqE; P: in Taylor; H: in out Fourier); -- H_err[I] := H_err[I]+P

  procedure AddBall(S: in Scalar; H: in out Fourier);                 -- Add |S| to H_err[0]_err
  procedure Copy(H1: in Fourier; H2: out Fourier);                    -- H2 := H1
  procedure ZeroConstant(H: in out Fourier);                          -- H_{0,0,0} := 0
  procedure QDegConvert(L,W: in Positive; H: in out Fourier);         -- Enlarge H by reducing QDegs
  function Coeff(M,N: Integer; K: Natural; H: Fourier) return Scalar; -- Coeff := H_{M,N,K}

  function NormFactor(R: Radius; M,N: Integer;
                      Update: in Boolean := False) return Scalar;     -- Norm of function E_{M,N}
  function Norm(R: Weights; H: Fourier) return Scalar;                -- Norm := ||H||_R
  function Norm(H: Fourier) return Scalar;                            -- Norm := ||H||_R  where R=Rho(H)

  procedure Neg(H1: in Fourier; H2: out Fourier);                     -- H2 := -H1
  function "-"(H: Fourier) return Fourier;                            -- Return -H
  procedure Mult(S: in Scalar; H: in out Fourier);                    -- H := S*H
  procedure Prod(S: in Scalar; H1: in Fourier; H2: out Fourier);      -- H2 := S*H1
  function "*"(S: Scalar; H: Fourier) return Fourier;                 -- Return S*H
  procedure Add(H1: in Fourier; H2: in out Fourier);                  -- H2 := H2+H1
  procedure Sum(H1,H2: in Fourier; H3: out Fourier);                  -- H3 := H1+H2
  function "+"(H1,H2: Fourier) return Fourier;                        -- Return H1+H2
  procedure Sub(H1: in Fourier; H2: in out Fourier);                  -- H2 := H2-H1
  procedure Diff(H1,H2: in Fourier; H3: out Fourier);                 -- H3 := H1-H2
  function "-"(H1,H2: Fourier) return Fourier;                        -- Return H1-H2
  procedure MultAdd(S: in Scalar; H1: in Fourier; H2: in out Fourier); -- H2 := H2+S*H1

  procedure Iminus(L: in Degrees; H: in out Fourier);                 -- H := I^minus H
  function Iminus(L: Degrees; H: Fourier) return Fourier;             -- Iminus := I^minus H
  procedure Iplus(L: in Degrees; H: in out Fourier);                  -- H := I^plus H
  function Iplus(L: Degrees; H: Fourier) return Fourier;              -- Iplus := I^minus H

  procedure MultAdd(S: in Scalar; M,N: in Integer; P1: in Taylor;
                    H2: in Fourier; H3: in out Fourier;
               D: in Degrees := Dmax; Monom1: in Boolean := False);   -- H3 := H3+S*Monom(M,N,P1)*H2
  procedure MultAdd(S: in Scalar; H1,H2: in Fourier;
                    H3: in out Fourier; D: in Degrees := Dmax);       -- H3 := H3+S*H1*H2
  procedure Mult(H1: in Fourier; H2: in out Fourier);                 -- H2 := H1*H2
  procedure Prod(R: in Weights; H1,H2: in Fourier; H3: out Fourier;
                 D: in Degrees := Dmax);                              -- H3 := H1*H2   (I^minus)
  procedure Prod(R: in Weights; S: in Scalar; H1,H2: in Fourier;
                 H3: out Fourier; D: in Degrees := Dmax);             -- H3 := S*H1*H2   (I^minus)
  function "*"(H1,H2: Fourier) return Fourier;                        -- Return H1*H
  procedure Powers(R: in Weights; H: in Fourier; Hp: in FouVecPtr);   -- Hp := Powers_of_H
  procedure Inv(R: in Weights; H1: in Fourier; H2: in out Fourier);   -- H2 := 1/H1, given approx inverse
  procedure Scale(S: in Scalar; H1: in Fourier; H2: out Fourier);     -- H2 := H1(.,.,S.)/S
  procedure Scale(S: in Scalar; H: in out Fourier);                   -- H := H(.,.,S.)/S
  function Scale(S: Scalar; H: Fourier) return Fourier;               -- Scale := H(.,.,S.)/S
  procedure ScaleRem1(Eps: in Scalar; H1: in Fourier;
                                      H2: out Fourier);               -- H2 := (1+Eps)*H1(.,./(1+Eps))-H1
  procedure ScaleRem2(Eps: in Scalar; H1: in Fourier;                 -- H2 := (1+Eps)*H1(.,./(1+Eps))-H
                                      H2: out Fourier);               --       -Eps*H1+Eps*Z*(D_3 H1)

  procedure OmegaDq(R2: in Radius; H1: in Fourier; H2: out Fourier);  -- H2 := (Omega*Grad) H1
  procedure Der3(R3: in Radius; H1: in Fourier; H2: out Fourier);     -- H2 := D_3 H1
  function Der3(R3: Radius; H: Fourier) return Fourier;               -- Der3 := D_3 H
  procedure Gradient(R: in Weights; H: in Fourier; Hdq,Hdp: out Fourier); -- OmegaDq and Der3
  function DerDer3(H: Fourier) return Fourier;                        -- DerDer3 := D_3 D_3 H
  procedure Comp3(H: in Fourier; Hp3: in FouVecPtr; Hr: out Fourier); -- Hr := H(.,.,H3), Hp3 are powers of H3
  procedure Comp3(H,H3: in Fourier; Hr: out Fourier);                 -- Hr := H(.,.,H3)
  procedure Comp3(H3: in Fourier; H: in out Fourier);                 -- H := H(.,.,H3)
  function Comp3(H,H3: Fourier) return Fourier;                       -- Comp3 := H(.,.,H3)
  procedure Comp3(H3: in Fourier; H: in FouVecPtr);                   -- Comp3 for each component
  procedure CosSin(H: in Fourier; Hw1,Hw2,Hc,Hs: out Fourier);        -- Hc := Cos(H), Hs := Sin(H)
  function Gamma(S,T,X: Scalar) return Scalar;                        -- Return Exp((1-S-T)*X)*Cosh(S*X)/Cosh(X)
  function GammaMax(S,T,X0: Scalar) return Scalar;                    -- Return Max_{X>=X0}Gamma(S,T,X)
  procedure IncCS(R: in Weights; Hc1,Hs1: in Fourier;
                  Hwc,Hws: out Fourier; Hc2,Hs2: in out Fourier);     -- Used by CompOmegaTmu0
  procedure UHatEtc(Hq: in Fourier; Rp,R: in Weights;
                    SigHat,UHat: out Scalar);                         -- Used by CompOmegaTmu0
  procedure CompOmegaTmu0(R: in Weights; H,Hq: in Fourier;            -- Hr := a*H(T(.+Hq*Omega),b.)
                         Hr: out Fourier; DoErr: in Boolean := True); --       a=Golden^4, b=-1/Golden^2
  procedure CompOmegaTmu0(R: in Weights; H: in FouVecPtr;
        Hq: in Fourier; Hr: in FouVecPtr; DoErr: in Boolean := True); -- CompOmegaTmu0 for each component

  procedure ZxDer3(H1: in Fourier; H2: out Fourier);                  -- H2 := Z D_3 H1  where Z(q,z)=z
  procedure DqDot(S1,S2: in Scalar; H: in out Fourier);               -- H := (S1*D_1+S2*D_2) H
  procedure Gig(Sig: in Rep; H: in out Fourier);                      -- H := ((Omega*Grad)/(omega*Grad))H
  procedure Gig(Sig: in Rep; H1: in Fourier; H2: out Fourier);        -- H2 := ((Omega*Grad)/(omega*Grad))H1
  function Gig(Sig: in Rep; H: Fourier) return Fourier;               -- Return ((Omega*Grad)/(omega*Grad))H
  procedure Dig(Kap: in Rep; H: in out Fourier);                      -- H := (D_3/(omega*Grad))H
  procedure Dig(Kap: in Rep; H1: in Fourier; H2: out Fourier);        -- H2 := (D_3/(omega*Grad))H1
  function Dig(Kap: in Rep; H: Fourier) return Fourier;               -- Return (D_3/(omega*Grad))H

----- procedures that are implemented only numerically

  procedure NumInv(Err: in Radius; H1: in Fourier; H2: out Fourier);  -- H2 := 1/H1
  function NumInv(Err: Radius; H: Fourier) return Fourier;            -- Return 1/H
  procedure NumInvComp3(Err: in Radius; H1,H2: in Fourier;
                        H3: out Fourier);                             -- H3 solves H2(.,.,H3)=H1
  procedure NumFixComp3(Err: in Radius; H: in Fourier;
                        Hz: out Fourier);                             -- Hz solves -H(.,.,.+Hz)=Hz
  function NumFixComp3(Err: Radius; H: Fourier) return Fourier;       -- Return solution Hz of -H(.,.,.+Hz)=Hz
  procedure NumInvAplus1(Err: in Radius; L: in Degrees;
                         Hdq,Hdp,F: in Fourier; Psi: out Fourier);    -- Used in RG_Ops, inverts A+1
  function NumInvAplus1(Err: Radius; L: Degrees;
                         Hdq,Hdp,F: Fourier) return Fourier;          -- Used in RG_Ops, inverts A+1

private

  Ierr: IerrType := IerrDef;

  type Norm_Factors is
    record
      R: Radius := Zero;
      C: ScalMat(FreqE1,FreqE2);
    end record;

  SavedFactors: Norm_Factors;

  procedure UpdateSavedFactors(R: in Radius);
  procedure Normalize(M,N: in Integer; M3,N3: out Integer; Flip: out Boolean);
  procedure GammaGuessX(Err: in Radius; S,T,X0: in Rep; XL,XR: out Rep);
  procedure EffDeg(M: in Integer; H: in Fourier; Dn: out Integer);
  procedure EffDeg(H: in Fourier; Dm: out Integer);

end Fouriers.Ops;