File : fouriers-ops.adb
with Messages, Reps.Ops;
use Messages;
pragma Elaborate_All (Messages, Reps.Ops);
package body Fouriers.Ops is
use Components_Ops;
Small: constant Rep := Half**10;
ScalOne: constant Scalar := Scal(One);
ScalTwo: constant Scalar := Scal(Two);
ScalHalf: constant Scalar := Scal(Half);
function IMin(I1,I2: Integer) return Integer renames Integer'Min;
function IMax(I1,I2: Integer) return Integer renames Integer'Max;
function Deg(P: Taylor) return Degree renames Components.Deg; -- want to inline it
pragma Inline (Imin,Imax,Deg);
function DmaxCount return Integer is
D: Integer := 0;
begin
for M in Freq1 loop
for N in Nmin(M) .. Nmax(M) loop
D := D+Dmax(M,N)+1;
end loop;
end loop;
return D;
end DmaxCount;
function ErrCount return Integer is
C: Integer;
begin
C := (Pdeg+1)*NumErrs(Lmax,Wmax);
for M in Freq1 loop
for N in Nmin(M) .. Nmax(M) loop
if Dmax(M,N) >= 0 then C := C+1; end if;
end loop;
end loop;
return C;
end ErrCount;
function InRange(M,N: Integer; K: Natural := 0; D: Degrees := Dmax) return Boolean is
begin
return M in Freq1 and then N in Freq2 and then K<=D(M,N);
end InRange;
procedure Component(K: Natural; H1: in Fourier; H2: out Fourier) is
E,S: Scalar;
begin
H2.P := H1.P;
H2.J := H1.J;
H2.R := H1.R;
SetZero(H2.C(0,0));
S := ScalZero;
for M in Freq1 loop
for N in Nmin(M,H1.J) .. Nmax(M) loop
ComponentSplit(K,H1.C(M,N),H2.C(M,N),E);
if not H1.P then
S := S+E*NormFactor(H1.R.Q,M,N,True);
end if;
end loop;
end loop;
if Trunc then return; end if;
if H1.P then
ResetErrComp(H2);
return;
end if;
for I in FreqE loop
ComponentSplit(K,H1.E(I),H2.E(I),E);
S := S+E;
end loop;
AddErrComp(Abs(S),H2.E(0));
end Component;
procedure SetComponent(M: in Freq1; N: in Freq2; K: in Natural; S: in Scalar; H: in out Fourier) is
--- allowing K's larger than Dmax(M,N)
begin
if InRange(M,N) and (H.J=0 or else M /= 0 or else N /= 0 or else S=ScalZero) then
SetComponent(K,S,H.C(M,N));
if K<Dho or else Err0(S) then return; end if;
H.P := False; -- S may depend on p
else
Message("Fouriers.Ops.SetComponent",Index_Error);
end if;
end SetComponent;
procedure AddComponent(M: in Freq1; N: in Freq2; K: in Natural; S: in Scalar; H: in out Fourier) is
begin
if InRange(M,N) and (H.J=0 or else M /= 0 or else N /= 0 or else S=ScalZero) then
AddComponent(K,S,H.C(M,N));
if K<Dho or else Err0(S) then return; end if;
H.P := False; -- S may depend on p
else
Message("Fouriers.Ops.AddComponent",Index_Error);
end if;
end AddComponent;
procedure SetErrComp(I: in FreqE; K: in Natural; S: in Scalar; H: in out Fourier) is
begin
if Trunc then return; end if;
if K <= Pdeg then
SetComponent(K,Ball(S),H.E(I));
if S /= ScalZero then H.P := False; end if;
else
Message("Fouriers.Ops.SetErrComp",Index_Error);
end if;
end SetErrComp;
procedure AddErrComp(I: in FreqE; P: in Taylor; H: in out Fourier) is
begin
if Trunc or else Deg(P)<0 then return; end if;
MultAdd(Ball(One),P,H.E(I));
H.P := False;
end AddErrComp;
procedure AddBall(S: in Scalar; H: in out Fourier) is
begin
if Trunc or else S=ScalZero then return; end if;
AddErrComp(Abs(S),H.E(0));
H.P := False;
end AddBall;
procedure Copy(H1: in Fourier; H2: out Fourier) is
begin
H2.P := H1.P;
H2.J := H1.J;
H2.R := H1.R;
for M in Freq1 loop
for N in Nmin(M) .. Nmax(M) loop
Copy(H1.C(M,N),H2.C(M,N));
end loop;
end loop;
if Trunc then return; end if;
for I in FreqE loop
Copy(H1.E(I),H2.E(I));
end loop;
end Copy;
procedure ZeroConstant(H: in out Fourier) is
begin
if Dho>0 then
SetComponent(0,0,0,ScalZero,H);
else
ResetCenter(H.C(0,0),0,0);
end if;
end ZeroConstant;
procedure QDegConvert(L,W: in Positive; H: in out Fourier) is
Lt: constant Positive := Imin(L,Lmax);
Wt: constant Positive := Imin(W,Wmax);
D: Integer;
begin
for M in Freq1 loop
for N in Nmin(M) .. Nmax(M) loop
Cleanup(H.C(M,N));
if Kcut(M,N,Lt,Wt)<0 and then Deg(H.C(M,N)) >= 0 then
if not Trunc then
MultAdd(Ball(NormFactor(H.R.Q,M,N)),H.C(M,N),H.E(Ierr(M,N)));
H.P := False;
end if;
SetZero(H.C(M,N));
end if;
end loop;
end loop;
if H.P then return; end if;
D := NumErrs(Lt,Wt);
for I in D+1 .. SizeE loop
Add(H.E(I),H.E(D));
SetZero(H.E(I));
end loop;
end QDegConvert;
function Coeff(M,N: Integer; K: Natural; H: Fourier) return Scalar is
C: Scalar := ScalZero;
begin
if H.P then
if InRange(M,N) then
C := Component(K,H.C(M,N));
end if;
else
for I in 0 .. Ierr(M,N) loop
C := C+Coeff(H.R.P,K,H.E(I));
end loop;
C := Ball(C/NormFactor(H.R.Q,M,N));
if InRange(M,N) then
C := C+Coeff(H.R.P,K,H.C(M,N));
end if;
end if;
return C;
end Coeff;
procedure UpdateSavedFactors(R: in Radius) is
Fm,Fn,Pm,Pmn: Scalar;
begin
if R=SavedFactors.R then return; end if;
SavedFactors.R := R;
Fn := Exp((-R)*InvGold);
Fm := Exp(Scal(R)); ---- start nonnegative m's
Pm := ScalOne;
for M in 0 .. SizeE1 loop
Pmn := Pm;
for N in FreqE2 loop
SavedFactors.C(M,N) := Half*(Pmn+Inv(Pmn));
Pmn := Pmn*Fn;
end loop;
Pm := Pm*Fm;
end loop;
Fm := Inv(Fm); ---- now negative m's
Pm := Fm;
for M in reverse -SizeE1 .. -1 loop
Pmn := Pm;
for N in FreqE2 loop
SavedFactors.C(M,N) := Half*(Pmn+Inv(Pmn));
Pmn := Pmn*Fn;
end loop;
Pm := Pm*Fm;
end loop;
end UpdateSavedFactors;
function NormFactor(R: Radius; M,N: Integer; Update: in Boolean := False) return Scalar is
begin
if N<0 then
return NormFactor(R,-M,-N,Update);
elsif R=SavedFactors.R then
if M in FreqE1 and then N in FreqE2 then
return SavedFactors.C(M,N);
end if;
elsif Update then
UpdateSavedFactors(R);
end if;
return Cosh(R*(Scal(M)+Rep(-N)*InvGold));
end NormFactor;
function Norm(R: Weights; H: Fourier) return Scalar is
E: Scalar := ScalZero;
begin
if not H.P then
for I in FreqE loop
E := E+Norm(R.P,H.E(I));
end loop;
end if;
if R.Q=SavedFactors.R then
for M in Freq1 loop
for N in Nmin(M,H.J) .. Nmax(M) loop
E := E+SavedFactors.C(M,N)*Norm(R.P,H.C(M,N));
end loop;
end loop;
return E;
end if;
declare
Fm,Fn,Pm,Pmn: Scalar;
begin
Fn := Exp((-R.Q)*InvGold);
Fm := Exp(Scal(R.Q)); ---- start nonnegative m's
Pm := ScalOne;
for M in 0 .. Size1 loop
Pmn := Pm;
for N in 0 .. Nmax(M) loop
E := E+Half*(Pmn+Inv(Pmn))*Norm(R.P,H.C(M,N));
Pmn := Pmn*Fn;
end loop;
Pm := Pm*Fm;
end loop;
Fm := Inv(Fm); ---- now negative m's
Pm := Fm;
for M in reverse -Size1 .. -1 loop
Pmn := Pm*Fn;
for N in 1 .. Nmax(M) loop
E := E+Half*(Pmn+Inv(Pmn))*Norm(R.P,H.C(M,N));
Pmn := Pmn*Fn;
end loop;
Pm := Pm*Fm;
end loop;
end;
return E;
end Norm;
function Norm(H: Fourier) return Scalar is
begin
return Norm(H.R,H);
end Norm;
procedure Neg(H1: in Fourier; H2: out Fourier) is
begin
H2.P := H1.P;
H2.J := H1.J;
H2.R := H1.R;
for M in Freq1 loop
for N in Nmin(M) .. Nmax(M) loop
Neg(H1.C(M,N),H2.C(M,N));
end loop;
end loop;
if Trunc then return; end if;
for I in FreqE loop
Copy(H1.E(I),H2.E(I));
end loop;
end Neg;
function "-"(H: Fourier) return Fourier is
Ht: Fourier;
begin
Neg(H,Ht);
return Ht;
end "-";
procedure Mult(S: in Scalar; H: in out Fourier) is
begin
if S = ScalZero then
SetZero(H.J,H.R,H);
return;
end if;
H.P := H.P and Err0(S); -- S may depend on p
for M in Freq1 loop
for N in Nmin(M,H.J) .. Nmax(M) loop
Mult(S,H.C(M,N));
end loop;
end loop;
if H.P then return; end if;
for I in FreqE loop
Mult(Abs(S),H.E(I));
end loop;
end Mult;
procedure Prod(S: in Scalar; H1: in Fourier; H2: out Fourier) is
begin
if S = ScalZero then
SetZero(H1.J,H1.R,H2);
return;
end if;
H2.P := H1.P and Err0(S); -- S may depend on p
H2.J := H1.J;
H2.R := H1.R;
for M in Freq1 loop
for N in Nmin(M) .. Nmax(M) loop
Prod(S,H1.C(M,N),H2.C(M,N));
end loop;
end loop;
if Trunc then return; end if;
if H1.P then
ResetErrComp(H2);
return;
end if;
for I in FreqE loop
Prod(Abs(S),H1.E(I),H2.E(I));
end loop;
end Prod;
function "*"(S: Scalar; H: Fourier) return Fourier is
Ht: Fourier;
begin
Prod(S,H,Ht);
return Ht;
end "*";
procedure Add(H1: in Fourier; H2: in out Fourier) is
begin
if H1.J /= H2.J then
Message("Fouriers.Ops.Add",Parity_Error);
end if;
H2.R := Min(H1.R,H2.R);
for M in Freq1 loop
for N in Nmin(M,H1.J) .. Nmax(M) loop
Add(H1.C(M,N),H2.C(M,N));
end loop;
end loop;
if H1.P then return; end if;
H2.P := Trunc;
for I in FreqE loop
Add(H1.E(I),H2.E(I));
end loop;
end Add;
procedure Sum(H1,H2: in Fourier; H3: out Fourier) is
begin
if H1.J /= H2.J then
Message("Fouriers.Ops.Sum",Parity_Error);
end if;
H3.P := H1.P and H2.P;
H3.J := H1.J;
H3.R := Min(H1.R,H2.R);
for M in Freq1 loop
for N in Nmin(M) .. Nmax(M) loop
Sum(H1.C(M,N),H2.C(M,N),H3.C(M,N));
end loop;
end loop;
if Trunc then return; end if;
if H3.P then
ResetErrComp(H3);
return;
end if;
for I in FreqE loop
Sum(H1.E(I),H2.E(I),H3.E(I));
end loop;
end Sum;
function "+"(H1,H2: Fourier) return Fourier is
Hr: Fourier;
begin
Sum(H1,H2,Hr);
return Hr;
end "+";
procedure Sub(H1: in Fourier; H2: in out Fourier) is
begin
if H1.J /= H2.J then
Message("Fouriers.Ops.Sub",Parity_Error);
end if;
H2.R := Min(H1.R,H2.R);
for M in Freq1 loop
for N in Nmin(M,H1.J) .. Nmax(M) loop
Sub(H1.C(M,N),H2.C(M,N));
end loop;
end loop;
if H1.P then return; end if;
H2.P := Trunc;
for I in FreqE loop
Add(H1.E(I),H2.E(I));
end loop;
end Sub;
procedure Diff(H1,H2: in Fourier; H3: out Fourier) is
begin
if H1.J /= H2.J then
Message("Fouriers.Ops.Diff",Parity_Error);
end if;
H3.P := H1.P and H2.P;
H3.J := H1.J;
H3.R := Min(H1.R,H2.R);
for M in Freq1 loop
for N in Nmin(M) .. Nmax(M) loop
Diff(H1.C(M,N),H2.C(M,N),H3.C(M,N));
end loop;
end loop;
if Trunc then return; end if;
if H3.P then
ResetErrComp(H3);
return;
end if;
for I in FreqE loop
Sum(H1.E(I),H2.E(I),H3.E(I));
end loop;
end Diff;
function "-"(H1,H2: Fourier) return Fourier is
Hr: Fourier;
begin
Diff(H1,H2,Hr);
return Hr;
end "-";
procedure MultAdd(S: in Scalar; H1: in Fourier; H2: in out Fourier) is
begin
if S=ScalZero then return; end if;
if H1.J /= H2.J then
Message("Fouriers.Ops.MultAdd",Parity_Error);
end if;
H2.P := Err0(S) and H1.P and H2.P; -- S may depend on p
H2.R := Min(H1.R,H2.R);
for M in Freq1 loop
for N in Nmin(M,H1.J) .. Nmax(M) loop
MultAdd(S,H1.C(M,N),H2.C(M,N));
end loop;
end loop;
if H1.P then return; end if;
for I in FreqE loop
MultAdd(Abs(S),H1.E(I),H2.E(I));
end loop;
end MultAdd;
procedure Iminus(L: in Degrees; H: in out Fourier) is
begin
for M in Freq1 loop
for N in Nmin(M,H.J) .. Nmax(M) loop
if L(M,N)<0 then
SetZero(H.C(M,N));
else
ResetComponent(H.C(M,N),L(M,N)+1);
end if;
end loop;
end loop;
end Iminus;
function Iminus(L: Degrees; H: Fourier) return Fourier is
Hr: Fourier;
begin
Copy(H,Hr);
Iminus(L,Hr);
return Hr;
end Iminus;
procedure Iplus(L: in Degrees; H: in out Fourier) is
begin
for M in Freq1 loop
for N in Nmin(M,1) .. Nmax(M) loop
if H.P then
ResetComponent(H.C(M,N),0,L(M,N));
else
ResetComponent(H.C(M,N),0,Imin(Dho-1,L(M,N)));
ResetCenter(H.C(M,N),Imin(Dho-1,L(M,N))+1,L(M,N));
end if;
end loop;
end loop;
end Iplus;
function Iplus(L: Degrees; H: Fourier) return Fourier is
Hr: Fourier;
begin
Copy(H,Hr);
Iplus(L,Hr);
return Hr;
end Iplus;
procedure Normalize(M,N: in Integer; M3,N3: out Integer; Flip: out Boolean) is
begin
if N<0 or else (N=0 and then M<0) then
M3 := -M; N3 := -N; Flip := True;
else
M3 := M; N3 := N; Flip := False;
end if;
end Normalize;
pragma Inline (Normalize);
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) is
J1: constant Parity := (H2.J+H3.J) mod 2;
Flip: Boolean;
M1,M3,N1,N3,I1: Integer;
Cp,Cm1,Cm2,F: Scalar;
begin
Normalize(M,N,M1,N1,Flip); --- (M,N) need not be inrange
if Flip and then J1=1 then
Cp := (-Half)*S;
else
Cp := Half*S;
end if;
Cm1 := Cp;
Cm2 := Cp;
if J1 < H2.J then
Cm1 := -Cp;
elsif J1 > H2.J then
Cm2 := -Cp;
elsif J1 = 1 then
Cp := -Cp;
end if;
if not Trunc then
if Min(H2.R,H3.R) /= H3.R then
Message("Fouriers.Ops.MultAdd",Domain_Violation);
end if;
H3.P := False;
UpdateSavedFactors(H3.R.Q);
end if;
for M2 in Freq1 loop
for N2 in Nmin(M2,H2.J) .. Nmax(M2) loop
if Deg(H2.C(M2,N2))>=0 then
Normalize(M1+M2,N1+N2,M3,N3,Flip);
if InRange(M3,N3) then --- frequency sum
MultAdd(H3.R.P,Cp,P1,H2.C(M2,N2),H3.C(M3,N3),D(M3,N3),Monom1);
elsif not Trunc then
F := Ball(Cp*SavedFactors.C(M3,N3));
MultAdd(H3.R.P,F,P1,H2.C(M2,N2),H3.E(Ierr(M3,N3)),Pdeg,Monom1);
end if;
Normalize(M1-M2,N1-N2,M3,N3,Flip);
if InRange(M3,N3) then --- frequency difference
if not Flip then --- 1st minus 2nd frequency
MultAdd(H3.R.P,Cm1,P1,H2.C(M2,N2),H3.C(M3,N3),D(M3,N3),Monom1);
else --- 2nd minus 1st frequency
MultAdd(H3.R.P,Cm2,P1,H2.C(M2,N2),H3.C(M3,N3),D(M3,N3),Monom1);
end if;
elsif not Trunc then
F := Ball(Cp*SavedFactors.C(M3,N3));
MultAdd(H3.R.P,F,P1,H2.C(M2,N2),H3.E(Ierr(M3,N3)),Pdeg,Monom1);
end if;
end if;
end loop;
end loop;
if H3.J=1 then
SetZero(H3.C(0,0));
end if;
if H2.P then return; end if;
if M1=0 and then N1=0 then I1 := -1; else I1 := Ierr(M1,N1); end if;
F := Ball(S*SavedFactors.C(M1,N1));
for I2 in FreqE loop
MultAdd(H3.R.P,F,P1,H2.E(I2),H3.E(Imax(0,I2-I1-1)),Pdeg,Monom1); --- add coeffs*errors
end loop;
end MultAdd;
procedure MultAdd(S: in Scalar; H1,H2: in Fourier; H3: in out Fourier; D: in Degrees := Dmax) is
I2: Integer;
F: Scalar;
begin
if not Trunc then
if Min(Min(H1.R,H2.R),H3.R) /= H3.R then
Message("Fouriers.Ops.MultAdd",Domain_Violation);
end if;
UpdateSavedFactors(H3.R.Q);
end if;
for M in Freq1 loop
for N in Nmin(M,H1.J) .. Nmax(M) loop
if Deg(H1.C(M,N))>=0 then
MultAdd(S,M,N,H1.C(M,N),H2,H3,D); --- includes coeffs*errors
end if;
end loop;
end loop;
if H1.P then return; end if;
UpdateSavedFactors(H3.R.Q);
for M in Freq1 loop
for N in Nmin(M,H2.J) .. Nmax(M) loop
if M=0 and then N=0 then I2 := -1; else I2 := Ierr(M,N); end if;
F := Ball(S*SavedFactors.C(M,N));
for I1 in FreqE loop
MultAdd(H3.R.P,F,H1.E(I1),H2.C(M,N),H3.E(Imax(0,I1-I2-1)),Pdeg); --- add errors*coeffs
end loop;
end loop;
end loop;
if H2.P then return; end if;
declare
E1,E2: Taylor;
begin
SetZero(E1);
SetZero(E2);
for I in FreqE loop
Add(H1.E(I),E1);
Add(H2.E(I),E2);
end loop;
MultAdd(H3.R.P,Abs(S),E1,E2,H3.E(0),Pdeg); --- add errors*errors
end;
end MultAdd;
procedure Mult(H1: in Fourier; H2: in out Fourier) is
Ht: Fourier;
begin
Copy(H2,Ht);
SetZero((H1.J+Ht.J) mod 2,Min(H1.R,Ht.R),H2);
MultAdd(ScalOne,H1,Ht,H2);
end Mult;
procedure Prod(R: in Weights; S: in Scalar; H1,H2: in Fourier; H3: out Fourier; D: in Degrees := Dmax) is
begin
SetZero((H1.J+H2.J) mod 2,R,H3);
MultAdd(S,H1,H2,H3,D);
end Prod;
procedure Prod(R: in Weights; H1,H2: in Fourier; H3: out Fourier; D: in Degrees := Dmax) is
begin
SetZero((H1.J+H2.J) mod 2,R,H3);
MultAdd(ScalOne,H1,H2,H3,D);
end Prod;
function "*"(H1,H2: Fourier) return Fourier is
Hr: Fourier;
begin
SetZero((H1.J+H2.J) mod 2,Min(H1.R,H2.R),Hr);
MultAdd(ScalOne,H1,H2,Hr);
return Hr;
end "*";
procedure Powers(R: in Weights; H: in Fourier; Hp: in FouVecPtr) is
I: Integer;
begin
if not (H.P or else Min(R,H.R)=R) then
Message("Fouriers.Ops.Powers",Domain_Violation);
end if;
Copy(H,Hp(1));
Hp(1).R := R;
I := 1;
for K in 2 .. Hp'Last loop
Prod(R,Hp(I),Hp(K-I),Hp(K));
if I+I=K then I := K; end if;
end loop;
end Powers;
procedure Inv(R: in Weights; H1: in Fourier; H2: in out Fourier) is
--- just adds ball around approximate inverse H2
begin
if Trunc then return; end if;
declare
Ea,Eb,Ec,Ed: Scalar;
Ht: Fourier;
begin
Prod(R,H1,H2,Ht);
AddComponent(0,0,0,-ScalOne,Ht);
Ea := Norm(Ht);
Eb := Norm(R,H1);
Ec := Norm(R,H2)*Ea;
Ed := Two*(Ea+Eb*Ec);
if Sup(Ed)>Half or Sup(Ed+Two*Eb*Ec)>=One then
Message("Fouriers.Ops.Inv",Domain_Violation); --- to stop, use Domain_Error
end if;
H2.P := False;
AddErrComp(Ec*Inv(ScalOne-Ed),H2.E(0));
end;
end Inv;
procedure Scale(S: in Scalar; H1: in Fourier; H2: out Fourier) is
begin
H2.P := H1.P; -- S is constant
H2.J := H1.J;
H2.R := H1.R;
H2.R.P := Inf(Scal(H2.R.P)/Abs(S));
for M in Freq1 loop
for N in Nmin(M) .. Nmax(M) loop
Scale(S,H1.C(M,N),H2.C(M,N));
end loop;
end loop;
if Trunc then return; end if;
if H1.P then
ResetErrComp(H2);
return;
end if;
for I in FreqE loop
Scale(Abs(S),H1.E(I),H2.E(I));
end loop;
end Scale;
procedure Scale(S: in Scalar; H: in out Fourier) is
begin
H.R.P := Inf(Scal(H.R.P)/Abs(S));
for M in Freq1 loop
for N in Nmin(M,H.J) .. Nmax(M) loop
Scale(S,H.C(M,N));
end loop;
end loop;
if H.P then return; end if;
for I in FreqE loop
Scale(Abs(S),H.E(I));
end loop;
end Scale;
function Scale(S: Scalar; H: Fourier) return Fourier is
Hr: Fourier;
begin
Scale(S,H,Hr);
return Hr;
end Scale;
procedure ScaleRem1(Eps: in Scalar; H1: in Fourier; H2: out Fourier) is
begin
SetZero(H1.J,H1.R,H2);
H2.P := H1.P; -- Eps is a constant
for M in Freq1 loop
for N in Nmin(M,H1.J) .. Nmax(M) loop
ScaleRem1(Eps,H1.C(M,N),H2.C(M,N),H1.P);
end loop;
end loop;
end ScaleRem1;
procedure ScaleRem2(Eps: in Scalar; H1: in Fourier; H2: out Fourier) is
begin
SetZero(H1.J,H1.R,H2);
H2.P := H1.P; -- Eps is a constant
for M in Freq1 loop
for N in Nmin(M,H1.J) .. Nmax(M) loop
ScaleRem2(Eps,H1.C(M,N),H2.C(M,N),H1.P);
end loop;
end loop;
end ScaleRem2;
procedure OmegaDq(R2: in Radius; H1: in Fourier; H2: out Fourier) is
S: Scalar;
begin
H2.P := H1.P;
H2.J := 1-H1.J;
H2.R.Q := R2;
H2.R.P := H1.R.P;
for M in Freq1 loop
for N in Nmin(M,1) .. Nmax(M) loop
if H1.J = 0 then
Prod(Scal(-M)+Rep(N)*InvGold,H1.C(M,N),H2.C(M,N));
else
Prod(Scal(M)+Rep(-N)*InvGold,H1.C(M,N),H2.C(M,N));
end if;
end loop;
end loop;
SetZero(H2.C(0,0));
if Trunc then return; end if;
if H1.P then
ResetErrComp(H2);
return;
end if; -- not really used beyond this
if R2>=H1.R.Q then
Message("Fouriers.Ops.OmegaDq",Domain_Violation);
end if;
S := GammaMax(Scal(R2)/H1.R.Q,ScalZero,ScalZero);
S := S/(Exp(ScalOne)*(Scal(H1.R.Q)-Scal(R2)));
for I in FreqE loop
Prod(Abs(S),H1.E(I),H2.E(I));
end loop;
end OmegaDq;
procedure Der3(R3: in Radius; H1: in Fourier; H2: out Fourier) is
Mt: Integer;
Ft: Scalar;
begin
H2.P := H1.P;
H2.J := H1.J;
H2.R := H1.R;
for M in Freq1 loop
for N in Nmin(M) .. Nmax(M) loop
Der(H1.R.P,R3,H1.C(M,N),H2.C(M,N),H1.P);
end loop;
end loop;
H2.R.P := R3;
if Trunc then return; end if;
if H1.P then
ResetErrComp(H2);
end if;
Fmax1(Scal(R3)/H1.R.P,Mt,Ft);
Ft := Ft/R3;
for I in FreqE loop
Prod(Ft,H1.E(I),H2.E(I));
end loop;
end Der3;
function Der3(R3: Radius; H: Fourier) return Fourier is
Hr: Fourier;
begin
Der3(R3,H,Hr);
return Hr;
end Der3;
procedure Gradient(R: in Weights; H: in Fourier; Hdq,Hdp: out Fourier) is
begin
OmegaDq(R.Q,H,Hdq);
Der3(R.P,H,Hdp);
end Gradient;
function DerDer3(H: Fourier) return Fourier is
Hr: Fourier;
begin
Copy(H,Hr);
for M in Freq1 loop
for N in Nmin(M) .. Nmax(M) loop
DerDer(Hr.C(M,N),Hr.P);
end loop;
end loop;
return Hr;
end DerDer3;
procedure Comp3(H: in Fourier; Hp3: in FouVecPtr; Hr: out Fourier) is
--- Hp3(K) is assumed to be the K-th power of Hp3(1)
R: Weights;
Hk: Fourier;
begin
TraceEnter("Comp3");
if Hp3(1).J /= 0 then
Message("Fouriers.Ops.Comp3",Not_Implemented);
end if;
R.Q := Radius'Min(H.R.Q,Hp3(1).R.Q);
R.P := Hp3(1).R.P;
SetZero(H.J,R,Hr);
if not H.P then
R.P := Sup(Norm(R,Hp3(1)));
if R.P>H.R.P then
Message("Fouriers.Ops.Comp3",Domain_Violation);
end if;
end if;
for K in reverse 1 .. Pdeg loop
Component(K,H,Hk);
MultAdd(ScalOne,Hk,Hp3(K),Hr);
end loop;
Component(0,H,Hk);
Add(Hk,Hr);
TraceLeave;
end Comp3;
procedure Comp3(H,H3: in Fourier; Hr: out Fourier) is
Hp3: FouVecPtr := new FouVec(1 .. Pdeg);
begin
Powers((Radius'Min(H.R.Q,H3.R.Q),H3.R.P),H3,Hp3);
Comp3(H,Hp3,Hr);
end Comp3;
procedure Comp3(H3: in Fourier; H: in out Fourier) is
Hr: Fourier;
Hp3: FouVecPtr := new FouVec(1 .. Pdeg);
begin
Powers((Radius'Min(H.R.Q,H3.R.Q),H3.R.P),H3,Hp3);
Comp3(H,Hp3,Hr);
Copy(Hr,H);
end Comp3;
function Comp3(H,H3: Fourier) return Fourier is
Hr: Fourier;
Hp3: FouVecPtr := new FouVec(1 .. Pdeg);
begin
Powers((Radius'Min(H.R.Q,H3.R.Q),H3.R.P),H3,Hp3);
Comp3(H,Hp3,Hr);
return Hr;
end Comp3;
procedure Comp3(H3: in Fourier; H: in FouVecPtr) is
Hr: Fourier;
Hp3: FouVecPtr := new FouVec(1 .. Pdeg);
begin
Powers((Radius'Min(H(1).R.Q,H3.R.Q),H3.R.P),H3,Hp3);
for L in H'Range loop
Comp3(H(L),Hp3,Hr);
Copy(Hr,H(L));
end loop;
end Comp3;
procedure CosSin(H: in Fourier; Hw1,Hw2,Hc,Hs: out Fourier) is
--- uses Hw1 and Hw2 as workspace only
Iter: constant Integer := 6;
Small: constant Rep := Half;
K: Natural := 0;
S,E,Ec,Es: Scalar;
begin
E := Norm(H);
while Sup(E)>Small loop
K := K+1;
E := Half*E;
end loop;
SetZero(0,H.R,Hc);
AddComponent(0,0,0,ScalOne,Hc);
S := Scal(Half**K);
Prod(S,H,Hs);
Copy(Hs,Hw1);
Ec := E**(2*Iter+1);
for I in 1 .. Iter loop -- Taylor polynomials
Prod(H.R,S,H,Hw1,Hw2);
Mult(Inv(Scal(-2*I)),Hw2);
Add(Hw2,Hc);
Prod(H.R,S,H,Hw2,Hw1);
Mult(Inv(Scal(2*I+1)),Hw1);
Add(Hw1,Hs);
Ec := Ec/Rep(2*I*(2*I+1));
end loop;
if not Trunc then
Ec := E*Ec/Rep(2*Iter+2);
Es := E*Ec/Rep(2*Iter+3);
Ec := Ec/(ScalOne-Sqr(E)/Rep((2*Iter+3)*(2*Iter+4)));
AddErrComp(Ec,Hc.E(0)); -- Remainder for Cos(H)
Es := Es/(ScalOne-Sqr(E)/Rep((2*Iter+4)*(2*Iter+5)));
AddErrComp(Es,Hs.E(0)); -- Remainder for Sin(H)
Hc.P := False;
Hs.P := False;
end if;
for I in 1 .. K loop -- double angle formula, K times
Prod(H.R,Hs,Hs,Hw1);
Prod(H.R,Hs,Hc,Hw2);
Prod(ScalTwo,Hw2,Hs);
Prod(H.R,Hc,Hc,Hw2); -- could avoid one multiplication
Diff(Hw2,Hw1,Hc);
end loop;
end CosSin;
procedure GammaGuessX(Err: in Radius; S,T,X0: in Rep; XL,XR: out Rep) is
--- Guess where Gamma takes its maximum value
use Reps.Ops;
Esx,Ex,XM,YL,YM,YR: Rep;
function DGammaSign(X: Rep) return Rep is
begin
Esx := Exp(-Two*S*X);
Ex:= Exp(-Two*X);
return Two*Ex*(One+Esx)/(One+Ex)-(T+Two*S)*Esx-T;
end DGammaSign;
begin
XR := X0;
YR := DGammaSign(XR);
if YR<=Zero then
XL := X0;
return;
end if;
while YR>Zero loop
XL := XR;
YL := YR;
XR := XL+One;
YR := DGammaSign(XR);
end loop;
while XR-XL>Err loop
XM := Half*(XL+XR);
YM := DGammaSign(XM);
if YM>Zero then
XL := XM;
YL := YM;
elsif YM<Zero then
XR := XM;
YR := YM;
else
XL := XM;
XR := XM;
return;
end if;
end loop;
end GammaGuessX;
function Gamma(S,T,X: Scalar) return Scalar is
begin
return Exp(-T*X)*(ScalOne+Exp((-Two)*S*X))/(ScalOne+Exp((-Two)*X));
end Gamma;
function GammaMax(S,T,X0: Scalar) return Scalar is
Err: constant Radius := 1.0E-6;
Eps: constant Radius := 1.0E-5;
XL,XR: Rep;
XM,Esx,Ex: Scalar;
function DGammaSign(X: Scalar) return Scalar is
begin
Esx := Exp((-Two)*S*X);
Ex:= Exp((-Two)*X);
return Two*Ex*(ScalOne+Esx)/(ScalOne+Ex)-(T+Two*S)*Esx-T;
end DGammaSign;
begin
if Inf(T)<Zero then
Message("RG_Ops.GammaMax",Parameter_Error);
end if;
if Sup(DGammaSign(X0))<Zero then
return Gamma(S,T,X0);
end if;
GammaGuessX(Err,Approx(S),Approx(T),Approx(X0),XL,XR);
XM := Scal(Half*(XR+XL))+Ball(Half*(XR-XL)+Eps);
if not Trunc and then (Inf(DGammaSign(Scal(Inf(XM))))<Zero or Sup(DGammaSign(Scal(Sup(XM))))>Zero) then
Message("RG_Ops.GammaMax",Parameter_Error);
end if;
return Gamma(S,T,XM);
end GammaMax;
procedure EffDeg(M: in Integer; H: in Fourier; Dn: out Integer) is
begin
Dn := -1;
for N in Nmin(M,H.J) .. Nmax(M) loop
if not IsZero(H.C(M,N)) then
Dn := Imax(N,Dn);
end if;
end loop;
end EffDeg;
procedure EffDeg(H: in Fourier; Dm: out Integer) is
begin
Dm := -1;
for M in Freq1 loop
for N in Nmin(M,H.J) .. Nmax(M) loop
if not IsZero(H.C(M,N)) then
Dm := Imax(Abs(M),Dm);
end if;
end loop;
end loop;
end EffDeg;
procedure IncCS(R: in Weights; Hc1,Hs1: in Fourier; Hwc,Hws: out Fourier; Hc2,Hs2: in out Fourier) is
--- make Hc2 := cos(H2+H1), Hs2 := sin(H2+H1)
--- from Hc1=cos(H1), Hs1=sin(H1), Hc2=cos(H2), Hs2=sin(H2)
--- uses Hwc and Hws as workspace only
begin
Copy(Hc2,Hwc);
Copy(Hs2,Hws);
Prod(R,Hc1,Hws,Hs2);
MultAdd(ScalOne,Hs1,Hwc,Hs2);
Prod(R,Hc1,Hwc,Hc2);
Multadd(-ScalOne,Hs1,Hws,Hc2);
end IncCS;
procedure UHatEtc(Hq: in Fourier; Rp,R: in Weights; SigHat,UHat: out Scalar) is
B: Weights;
begin
B.Q := Sup(Norm(R,Hq));
B.P := R.P;
SigHat := Scal(Rp.Q)-R.Q*InvGold-B.Q*InvGold;
if Inf(Rp.P*Sqr(Gold))<B.P then
Message("Fouriers.Ops.UHatEtc",Domain_Violation);
end if;
UHat := GammaMax(R.Q*InvGold/Rp.Q,SigHat/Rp.Q,ScalZero);
end UHatEtc;
procedure CompOmegaTmu0(R: in Weights; H,Hq: in Fourier; Hr: out Fourier; DoErr: in Boolean := True) is
S0: constant Scalar := -Gold-ScalOne; -- same as -Gold^2
S2: constant Scalar := -Sqr(InvGold); -- same as -1/Gold^2
Dm,Dn: Integer;
M: Freq1;
S1,SigHat,UHat: Scalar;
P: Taylor;
Hw1,Hw2,Hc1,Hs1,Hc2,Hs2,Hmc,Hms,Hmnc,Hmns: Fourier;
begin
TraceEnter("CompOmegaTmu0");
if Hq.J /= 1 then
Message("Fouriers.Ops.CompOmegaTmu0",Not_Implemented);
end if;
SetZero(H.J,R,Hr);
if DoErr and not H.P then -- could be improved with more info on H
for I in FreqE loop
Add(H.E(I),Hr.E(0));
end loop;
UHatEtc(Hq,H.R,R,SigHat,UHat); -- domains get checked here
Scale(Abs(S2),Hr.E(0));
Mult(UHat*Abs(S0),Hr.E(0));
end if;
EffDeg(H,Dm); -- looks only at H.C
if Dm<0 then
TraceLeave;
return;
end if;
S1 := Rep(2*H.J-1)*S0;
CosSin((-InvGold)*Hq,Hw1,Hw2,Hc1,Hs1);
CosSin(Sqr(InvGold)*Hq,Hw1,Hw2,Hc2,Hs2);
SetZero(0,R,Hmc);
SetComponent(0,0,0,ScalOne,Hmc);
SetZero(1,R,Hms);
for Mabs in 0 .. Dm loop
if Mabs>0 then
IncCS(R,Hc1,Hs1,Hw1,Hw2,Hmc,Hms);
end if;
for I in 0 .. 1 loop
M := (2*I-1)*Mabs;
Hms := -Hms;
if M<0 or else I=1 then
Copy(Hmc,Hmnc);
Copy(Hms,Hmns);
EffDeg(M,H,Dn);
for N in Nmin(M,H.J) .. Dn loop
if N > 0 then
IncCS(R,Hc2,Hs2,Hw1,Hw2,Hmnc,Hmns);
end if;
Scale(S2,H.C(M,N),P);
MultAdd(S0,N,M+N,P,Hmnc,Hr);
MultAdd(S1,N,M+N,P,Hmns,Hr);
end loop;
end if;
end loop;
end loop;
TraceLeave;
end CompOmegaTmu0;
procedure CompOmegaTmu0(R: in Weights; H: in FouVecPtr; Hq: in Fourier; Hr: in FouVecPtr; DoErr: in Boolean := True) is
--- R.P does not enter directly
Lf: constant Positive := Imin(H'Last,Hr'Last);
S0: constant Scalar := -Gold-ScalOne; -- same as -Gold^2
S2: constant Scalar := Gold-ScalTwo; -- same as -1/Gold^2
Dt,Dm,Dn: Integer;
M: Freq1;
S1,SigHat,UHat: Scalar;
P: Taylor;
Hw1,Hw2,Hc1,Hs1,Hc2,Hs2,Hmc,Hms,Hmnc,Hmns: Fourier;
begin
TraceEnter("CompOmegaTmu0");
if Hq.J /= 1 then
Message("Fouriers.Ops.CompOmegaTmu0",Not_Implemented);
end if;
for L in 1 .. Lf loop
SetZero(H(L).J,R,Hr(L));
if DoErr and not H(L).P then -- can be improved with more info on H
for I in FreqE loop
Add(H(L).E(I),Hr(L).E(0));
end loop;
UHatEtc(Hq,H(L).R,R,SigHat,UHat);
Scale(Abs(S2),Hr(L).E(0));
Mult(UHat*Abs(S0),Hr(L).E(0));
end if;
end loop;
Dm := -1;
for L in 1 .. Lf loop
EffDeg(H(L),Dt); -- looks only at H(L).C
Dm := IMax(Dm,Dt);
end loop;
if Dm<0 then
TraceLeave;
return;
end if;
CosSin((-InvGold)*Hq,Hw1,Hw2,Hc1,Hs1);
CosSin(Sqr(InvGold)*Hq,Hw1,Hw2,Hc2,Hs2);
SetZero(0,R,Hmc);
SetComponent(0,0,0,ScalOne,Hmc);
SetZero(1,R,Hms);
for Mabs in 0 .. Dm loop
if Mabs>0 then
IncCS(R,Hc1,Hs1,Hw1,Hw2,Hmc,Hms);
end if;
for I in 0 .. 1 loop
M := (2*I-1)*Mabs;
Hms := -Hms;
if M<0 or else I=1 then
Copy(Hmc,Hmnc);
Copy(Hms,Hmns);
Dn := -1;
for L in 1 .. Lf loop
EffDeg(M,H(L),Dt);
Dn := IMax(Dn,Dt);
end loop;
for N in Nmin(M) .. Dn loop
if N > 0 then
IncCS(R,Hc2,Hs2,Hw1,Hw2,Hmnc,Hmns);
end if;
for L in 1 .. Lf loop
Scale(S2,H(L).C(M,N),P);
S1 := Scal(2*H(L).J-1)*S0;
MultAdd(S0,N,M+N,P,Hmnc,Hr(L));
MultAdd(S1,N,M+N,P,Hmns,Hr(L));
end loop;
end loop;
end if;
end loop;
end loop;
TraceLeave;
end CompOmegaTmu0;
procedure ZxDer3(H1: in Fourier; H2: out Fourier) is
begin
H2.P := H1.P;
H2.J := H1.J;
H2.R := H1.R;
for M in Freq1 loop
for N in Nmin(M) .. Nmax(M) loop
ZxDer(H1.C(M,N),H2.C(M,N),H1.P);
end loop;
end loop;
if Trunc then return; end if;
ResetErrComp(H2);
end ZxDer3;
procedure DqDot(S1,S2: in Scalar; H: in out Fourier) is
begin
if not H.P then
Message("Fouriers.Ops.DqDot",Not_Implemented);
end if;
H.J := 1-H.J;
for M in Freq1 loop
for N in Nmin(M,1) .. Nmax(M) loop
if H.J = 0 then
Mult(Rep(-M)*S1+Rep(-N)*S2,H.C(M,N));
else
Mult(Rep(M)*S1+Rep(N)*S2,H.C(M,N));
end if;
end loop;
end loop;
SetZero(H.C(0,0));
end DqDot;
procedure Gig(Sig: in Rep; H: in out Fourier) is
-- Sig defines a cutoff a la Iminus, unless H.P
begin
if not IsZero(H.C(0,0)) then
Message("Fouriers.Ops.Gig",Domain_Violation);
end if;
for M in Freq1 loop
for N in Nmin(M,1) .. Nmax(M) loop
if Deg(H.C(M,N))>=0 then
Mult((Scal(M)+Rep(-N)*InvGold)/(Rep(M)*InvGold+Scal(N)),H.C(M,N));
end if;
end loop;
end loop;
SetZero(H.C(0,0));
if H.P then return; end if;
for I in FreqE loop
Mult(ScalOne/Sig,H.E(I)); -- operator_norm <= 1/Sig
end loop;
end Gig;
procedure Gig(Sig: in Rep; H1: in Fourier; H2: out Fourier) is
-- Sig defines a cutoff a la Iminus, unless H1.P
begin
if not IsZero(H1.C(0,0)) then
Message("Fouriers.Ops.Gig",Domain_Violation);
end if;
H2.P := H1.P;
H2.J := H1.J;
H2.R := H1.R;
for M in Freq1 loop
for N in Nmin(M,1) .. Nmax(M) loop
if Deg(H1.C(M,N))>=0 then
Prod((Scal(M)+Rep(-N)*InvGold)/(Rep(M)*InvGold+Scal(N)),H1.C(M,N),H2.C(M,N));
else
SetZero(H2.C(M,N));
end if;
end loop;
end loop;
SetZero(H2.C(0,0));
if Trunc then return; end if;
if H1.P then
ResetErrComp(H2);
return;
end if;
for I in FreqE loop
Prod(ScalOne/Sig,H1.E(I),H2.E(I)); -- operator_norm <= 1/Sig
end loop;
end Gig;
function Gig(Sig: in Rep; H: Fourier) return Fourier is
Hr: Fourier;
begin
Gig(Sig,H,Hr);
return Hr;
end Gig;
procedure Dig(Kap: in Rep; H: in out Fourier) is
-- Kap defines a cutoff a la Iminus, unless H.P
KapInv: constant Radius := SupAbs(ScalOne/Kap);
begin
if not IsZero(H.C(0,0)) then
Message("Fouriers.Ops.Dig",Domain_Violation);
end if;
if H.J /= 0 then
Message("Fouriers.Ops.Dig",Not_Implemented);
end if;
H.J := 1;
for M in Freq1 loop
for N in Nmin(M,1) .. Nmax(M) loop
DerMult(H.R.P,KapInv,Inv(Rep(M)*InvGold+Scal(N)),H.C(M,N),H.P);
end loop;
end loop;
SetZero(H.C(0,0));
if H.P then return; end if;
for I in FreqE loop
DerMult(H.R.P,KapInv,ScalOne/Kap,H.E(I));
end loop;
end Dig;
procedure Dig(Kap: in Rep; H1: in Fourier; H2: out Fourier) is
-- Kap defines a cutoff a la Iminus, unless H1.P
KapInv: constant Radius := SupAbs(ScalOne/Kap);
begin
if not IsZero(H1.C(0,0)) then
Message("Fouriers.Ops.Dig",Domain_Violation);
end if;
if H1.J /= 0 then
Message("Fouriers.Ops.Dig",Not_Implemented);
end if;
H2.P := H1.P;
H2.J := 1;
H2.R := H1.R;
for M in Freq1 loop
for N in Nmin(M,1) .. Nmax(M) loop
DerProd(H1.R.P,KapInv,Inv(Rep(M)*InvGold+Scal(N)),H1.C(M,N),H2.C(M,N),H1.P);
end loop;
end loop;
SetZero(H2.C(0,0));
if Trunc then return; end if;
if H1.P then
ResetErrComp(H2);
return;
end if;
for I in FreqE loop
DerProd(H1.R.P,KapInv,ScalOne/Kap,H1.E(I),H2.E(I));
end loop;
end Dig;
function Dig(Kap: in Rep; H: Fourier) return Fourier is
Hr: Fourier;
begin
Dig(Kap,H,Hr);
return Hr;
end Dig;
----- procedures that are implemented only numerically
procedure NumInv(Err: in Radius; H1: in Fourier; H2: out Fourier) is
--- iterate h2 -> h2*[2-h1*h2]
S,E: Scalar;
Hw1,Hw2: Fourier;
begin
if not Trunc then
Message("Fouriers.Ops.NumInv",Numeric_Only);
end if;
S := Component(0,0,0,H1);
if S /= ScalZero then
Prod(-Inv(Sqr(S)),H1,H2);
AddComponent(0,0,0,ScalTwo/S,H2);
for I in 1 .. 30 loop
Prod(H1.R,H1,H2,Hw1);
AddComponent(0,0,0,-ScalOne,Hw1);
E := Norm(H1.R,Hw1);
AddComponent(0,0,0,-ScalOne,Hw1);
Prod(H1.R,H2,Hw1,Hw2);
Neg(Hw2,H2);
if Sup(E) <= Err then
return;
end if;
end loop;
Message("Fouriers.Ops.NumInv",Giving_Up);
else
Message("Fouriers.Ops.NumInv",Division_By_Zero);
end if;
end NumInv;
function NumInv(Err: Radius; H: Fourier) return Fourier is
Hr: Fourier;
begin
NumInv(Err,H,Hr);
return Hr;
end NumInv;
procedure NumInvComp3(Err: in Radius; H1,H2: in Fourier; H3: out Fourier) is
--- iterate h3 -> h3-[h2(..,h3)-h1]/(D_3 h2)(..,h3)
R: Weights;
E: Scalar;
H,Hd,Hd2,Hd3: Fourier;
begin
TraceEnter("NumInvComp3");
if not Trunc then
Message("Fouriers.Ops.NumInvComp3",Numeric_Only);
end if;
if H1.J /= H2.J then
Message("Fouriers.Ops.NumInvComp3",Not_Implemented);
end if;
R := Min(H1.R,H2.R); -- arbitrary
SetZero(0,R,H3);
InvComp(R.P,Err,H1.C(0,0),H2.C(0,0),H3.C(0,0));
Der3(R.P,H2,H);
NumInv(Err,H,Hd2); -- need not be that precise
for I in 1 .. 20 loop
Comp3(H2,H3,H);
Sub(H1,H);
Comp3(Hd2,H3,Hd);
Prod(R,H,Hd,Hd3);
Sub(Hd3,H3);
E := Norm(R,Hd3);
if Sup(E) <= Err then
TraceLeave;
return;
end if;
end loop;
Message("Fouriers.Ops.NumInvComp3",Giving_Up);
TraceLeave;
end NumInvComp3;
procedure NumFixComp3(Err: in Radius; H: in Fourier; Hz: out Fourier) is
H1,H2: Fourier;
begin
Copy(H,H2);
AddComponent(0,0,1,ScalOne,H2);
SetZero(0,H.R,H1);
SetComponent(0,0,1,ScalOne,H1);
NumInvComp3(Err,H1,H2,Hz);
-- AddComponent(0,0,1,-ScalOne,Hz); -- originally last line
Comp3(H,Hz,H1); -- only for contractions
Neg(H1,Hz); -- only for contractions
end NumFixComp3;
function NumFixComp3(Err: Radius; H: Fourier) return Fourier is
Hz: Fourier;
begin
NumFixComp3(Err,H,Hz);
return Hz;
end NumFixComp3;
procedure NumInvAplus1(Err: in Radius; L: in Degrees; Hdq,Hdp,F: in Fourier; Psi: out Fourier) is
--- solve (I-)[A+1]Psi=(I-)F by iteration Psi -> (I-)[F-A(Psi)]
--- where A(Psi)=Hdp*Gig(Psi)-Hdq*Dig(Psi)
Iter: constant Integer := 50;
I: Integer := 1;
E: Scalar;
D: Degrees(Freq1,Freq2);
Psid,Psit: Fourier;
begin
if not Trunc then
Message("Fouriers.Ops.NumInvAplus1",Numeric_Only);
end if;
for M in Freq1 loop
for N in Nmin(M) .. Nmax(M) loop
D(M,N) := Imin(L(M,N),Dmax(M,N));
end loop;
end loop;
Copy(F,Psit);
Iminus(L,Psit); -- start with Psit := (I-)F
loop
Copy(F,Psi);
Dig(Small,Psit,Psid);
MultAdd(ScalOne,Hdq,Psid,Psi,D); -- add (I-)*Hdq*Dig(Psit)
Gig(Small,Psit,Psid);
MultAdd(-ScalOne,Hdp,Psid,Psi,D); -- subtract (I-)*Hdp*Gig(Psit)
E := Norm(Psi-Psit);
exit when I=Iter or Sup(E)<Err;
Copy(Psi,Psit);
I := I+1;
end loop;
if I=Iter then
Message("RG_Ops.NumInvAplus1",Giving_Up);
end if;
Show("NumInvAplus1: ",I,False);
Show(" |Psi-Psit| =",Sup(E));
end NumInvAplus1;
function NumInvAplus1(Err: Radius; L: Degrees; Hdq,Hdp,F: Fourier) return Fourier is
Psi: Fourier;
begin
NumInvAplus1(Err,L,Hdq,Hdp,F,Psi);
return Psi;
end NumInvAplus1;
end Fouriers.Ops;