
package body FltChains is

  function Component(K: Integer; V: FltVec) return Flt is
    --- assuming V'First is 0 or -V'Last
    AbsK: constant Integer := abs(K);
  begin
    if AbsK>V'Last then return Zero; end if;
    if V'First<0 then return V(K); end if;
    if K<0 and then Par=1 then return -V(AbsK); end if;
    return V(AbsK);
  end Component;

  procedure Normalize(V: in out FltVec) is
    --- divide by max absolute value
    R: Radius := Zero;
  begin
    for K in V'Range loop R := RMax(abs(V(K)),R); end loop;
    for K in V'Range loop V(K) := V(K)/R; end loop;
  end Normalize;

  function Translate(V: FltVec; T: Integer) return FltVec is
    use Flt_Vectors;
    U: FltVec(V'Range);
  begin
    SetZero(U);
    for K in U'Range loop
      U(K) := Component(K-T,V);
    end loop;
    return U;
  end Translate;

  procedure Symmetrize(V: in out FltVec; P: in Parity := Par) is
    --- assuming V'First is 0 or -V'Last
  begin
    if P=1 then V(0) := Zero; end if;
    if V'First<0 then
      if P=0 then
        for K in 1 .. V'Last loop V(K) := Half*(V(K)+V(-K)); V(-K) := V(K); end loop;
      else
        for K in 1 .. V'Last loop V(K) := Half*(V(K)-V(-K)); V(-K) := -V(K); end loop;
      end if;
    end if;
  end Symmetrize;

  procedure SumR(V: in FltChain; U: out FltVec) is
    -- solve U(t+1/2)+U(t-1/2)=V(t) starting from +infinity
    use Flt_Vectors;
    F: constant Integer := U'First;
    K: Integer;
    R: Flt;
  begin
    SetZero(U);
    for L in reverse U'Range loop
      R := Component(L,V);
      K := L-N;
      while K >= F loop
        U(K) := U(K)-R;
        K := K-2*N;
      end loop;
    end loop;
  end SumR;

  procedure TentConv(Width: in Positive; V1: in FltChain; V2: out FltChain) is
    --- convolution with tent function
    R: constant Radius := One/Flt(Width*Width);
  begin
    for K in V2'Range loop
      V2(K) := Zero;
      for M in 1-Width .. Width-1 loop
        V2(K) := V2(K)+Flt(Width-abs(M))*R*Component(K+M,V1);
      end loop;
    end loop;
    if Par=1 then V2(0) := Zero; end if;
  end TentConv;

  function Simple_Random(Decay: Radius) return FltChain is
    use Flt_EF;
    Width: constant Natural := N/2;
    Dec: constant Radius := Exp(Log(Decay)/Flt(2*N));
    Fac: Radius := One;
    V,W: FltChain;
  begin
    for K in V'Range loop
      V(K) := Fac*Simple_Random;
      Fac := Fac/Dec;
    end loop;
    if Par=1 then V(0) := Zero; end if;
    TentConv(Width,V,W);
    Normalize(W);
    return W;
  end Simple_Random;

  procedure PutPlot(F: in File_Type; V: in FltVec) is
    use FLt_IO;
    Digs: constant Positive := 4;
    L: constant Positive := V'Last;
    X,Y: Flt;
  begin
    for K in -L .. L loop
      X := Flt(K)/Flt(2*N);
      Put(F,XFac*X,3,Digs,4);
      Y := Component(K,V);
      Put(F,YFac*Y,3,Digs,4);
      New_Line(F);
    end loop;
  end PutPlot;

  procedure AA(Scale: in Positive; V1: in FltChain; V2: out FltChain) is
  begin
    TentConv(Scale*2*N,V1,V2);
  end AA;

end FltChains;
