File : taylors2.adb

with Ints, Ints.IO, Reps.IO, Messages;
use Ints, Ints.IO, Reps.IO, Messages;

pragma Elaborate (Ints);

package body Taylors2 is

  function InDomain(R: RadPair; F: Taylor2) return Boolean is
  begin
    return ((R(1) <= F.R(1)) and (R(2) <= F.R(2))) or else IsNumbers(F);
  end InDomain;

  procedure SetRho(R: in RadPair; F: in out Taylor2) is
  begin
    if not InDomain(R,F) then Error("Taylors2.SetRho domain error"); end if;
    F.R := R;
  end SetRho;

  function IsZero(F: Taylor2) return Boolean is
  begin
    return F.C(1,1)=ScalZero and then (F.C=Tay_Zero.C);
  end IsZero;

  function Zero_Taylor2 return Taylor2 is
    F: Taylor2;
  begin
    F.R := Some_RadPair;
    for I in Power loop
      for J in Power loop
        F.C(I,J) := ScalZero;
      end loop;
    end loop;
    return F;
  end Zero_Taylor2;

  procedure SetZero(R: in RadPair; F: out Taylor2) is
  begin
    F := Tay_Zero;
    F.R := R;
  end SetZero;

  procedure ZapConst(F: in out Taylor2) is
  begin
    ModNumbers(F.C(0,0));
  end ZapConst;

  procedure AddConst(C: in Scalar; F: in out Taylor2) is
  begin
    Add(C,F.C(0,0));
  end AddConst;

  procedure UnitVector(I,J: in Natural; R: in RadPair; F: out Taylor2) is
--- norm one
    S: constant Scalar := (Scal(R(1))**I)*(Scal(R(2))**J);
  begin
    SetZero(R,F);
    Add(Inv(S),F.C(I,J));
  end UnitVector;

  function QuasiDeg1(D: Integer; F: Taylor2) return Integer is
  begin
    for I in reverse 0 .. D loop
      for J in reverse 0 .. D-I loop
        if F.C(I,J) /= ScalZero then return I; end if;
      end loop;
    end loop;
    return -1;
  end QuasiDeg1;

  function QuasiDeg(F: Taylor2) return Integer is
  begin
    for D in reverse Power loop
      for I in reverse 0 .. D loop
        if F.C(I,D-I) /= ScalZero then return D; end if;
      end loop;
    end loop;
    return -1;
  end QuasiDeg;

  function IsNumbers(F: Taylor2) return Boolean is
  begin
    if not Trunc then
      for I in Power loop
        for J in 0 .. PDeg-I loop
          if not IsNumbers(F.C(I,J)) then return False; end if;
        end loop;
      end loop;
    end if;
    return True;
  end IsNumbers;

  function CenterDisk(F: Taylor2) return Taylor2 is
    G: Taylor2 := F;
  begin
    if not Trunc then
      for I in Power loop
        for J in 0 .. PDeg-I loop
          G.C(I,J) := CenterDisk(G.C(I,J));
        end loop;
      end loop;
    end if;
    return G;
  end CenterDisk;

  procedure ModNumbers(F: in out Taylor2) is
  begin
    if Trunc then
      F := Tay_Zero;
    else
      for I in Power loop
        for J in 0 .. PDeg-I loop
          ModNumbers(F.C(I,J));
        end loop;
      end loop;
    end if;
  end ModNumbers;

  procedure ModNumbers(F: in Taylor2; G: out Taylor2) is
  begin
    if Trunc then
      SetZero(F.R,G);
    else
      G := F;
      for I in Power loop
        for J in 0 .. PDeg-I loop
          ModNumbers(G.C(I,J));
        end loop;
      end loop;
    end if;
  end ModNumbers;

  function ModNumbers(F: Taylor2) return Taylor2 is
    G: Taylor2;
  begin
    ModNumbers(F,G);
    return G;
  end ModNumbers;

  function Norm(R: RadPair; F: Taylor2) return Scalar is
    S1,S2: Scalar;
  begin
    S1 := ScalZero;
    if not InDomain(R,F) then Error("Taylors2.Norm domain error"); end if;
    for I in reverse Power loop
      S2 := ScalZero;
      for J in reverse 0 .. PDeg-I loop
        S2 := R(2)*S2+Norm(F.C(I,J));
      end loop;
      S1 := R(1)*S1+S2;
    end loop;
    return Disk(S1);
  end Norm;

  function Norm(F: Taylor2) return Scalar is
  begin
    return Norm(F.R,F);
  end Norm;

  procedure ToBalls(F: in out Taylor2) is
  begin
    if Trunc then
      F.C := Tay_Zero.C;
    else
      for I in Power loop
        for J in 0 .. PDeg-I loop
          ToBalls(F.C(I,J));
        end loop;
      end loop;
    end if;
  end ToBalls;

  procedure Neg(F: in out Taylor2) is
  begin
    for I in Power loop
      for J in 0 .. PDeg-I loop
        Neg(F.C(I,J));
      end loop;
    end loop;
  end Neg;

  procedure Neg(F: in Taylor2; G: out Taylor2) is
  begin
    G := F;
    Neg(G);
  end Neg;

  function "-"(F: Taylor2) return Taylor2 is
    G: Taylor2 := F;
  begin
    Neg(G);
    return G;
  end "-";

  procedure Add(F: in Taylor2; G: in out Taylor2) is
  begin
    G.R := Min(F.R,G.R);
    for I in Power loop
      for J in 0 .. PDeg-I loop
        Add(F.C(I,J),G.C(I,J));
      end loop;
    end loop;
  end Add;

  procedure Sum(F,G: in Taylor2; H: out Taylor2) is
  begin
    H := F;
    Add(G,H);
  end Sum;

  function "+"(F,G: Taylor2) return Taylor2 is
    H: Taylor2 := F;
  begin
    Add(G,H);
    return H;
  end "+";

  procedure Sub(F: in Taylor2; G: in out Taylor2) is
  begin
    G.R := Min(F.R,G.R);
    for I in Power loop
      for J in 0 .. PDeg-I loop
        Sub(F.C(I,J),G.C(I,J));
      end loop;
    end loop;
  end Sub;

  procedure Diff(F,G: in Taylor2; H: out Taylor2) is
  begin
    H := F;
    Sub(G,H);
  end Diff;

  function "-"(F,G: Taylor2) return Taylor2 is
    H: Taylor2 := F;
  begin
    Sub(G,H);
    return H;
  end "-";

  procedure Mult(K: in Integer; F: in out Taylor2) is
  begin
    if K=0 then
      F.C := Tay_Zero.C;
    elsif K /= 1 then
      for I in Power loop
        for J in 0 .. PDeg-I loop
           Mult(K,F.C(I,J));
        end loop;
      end loop;
    end if;
  end Mult;

  procedure Prod(K: in Integer; F: in Taylor2; G: out Taylor2) is
  begin
    G := F;
    Mult(K,G);
  end Prod;

  function "*"(K: Integer; F: Taylor2) return Taylor2 is
    G: Taylor2 := F;
  begin
    Mult(K,G);
    return G;
  end "*";

  procedure MultAdd(K: in Integer; F: in Taylor2; G: in out Taylor2) is
  begin
    if K=1 then
      Add(F,G);
    elsif K=-1 then
      Sub(F,G);
    elsif K /= 0 then
      G.R := Min(F.R,G.R);
      for I in Power loop
        for J in 0 .. PDeg-I loop
          MultAdd(K,F.C(I,J),G.C(I,J));
        end loop;
      end loop;
    end if;
  end MultAdd;

  procedure Mult(S: in Scalar; F: in out Taylor2) is
  begin
    if S=ScalZero then
      F.C := Tay_Zero.C;
    else
      for I in Power loop
        for J in 0 .. PDeg-I loop
          Mult(S,F.C(I,J));
        end loop;
      end loop;
    end if;
  end Mult;

  procedure Prod(S: in Scalar; F: in Taylor2; G: out Taylor2) is
  begin
    G := F;
    Mult(S,G);
  end Prod;

  function "*"(S: Scalar; F: Taylor2) return Taylor2 is
    G: Taylor2 := F;
  begin
    Mult(S,G);
    return G;
  end "*";

  procedure MultAdd(S: in Scalar; F: in Taylor2; G: in out Taylor2) is
  begin
    if S /= ScalZero then
      G.R := Min(F.R,G.R);
      for I in Power loop
        for J in 0 .. PDeg-I loop
          MultAdd(S,F.C(I,J),G.C(I,J));
        end loop;
      end loop;
    end if;
  end MultAdd;

  procedure Mult(F: in Taylor2; G: in out Taylor2) is
    H: constant Taylor2 := G;
  begin
    G := Tay_Zero;
    MultAdd(F,H,G);
  end Mult;

  procedure Prod(F,G: in Taylor2; H: out Taylor2) is
  begin
    SetZero(F.R,H);
    MultAdd(F,G,H);
  end Prod;

  function "*"(F,G: Taylor2) return Taylor2 is
    H: Taylor2 := Tay_Zero;
  begin
    MultAdd(F,G,H);
    return H;
  end "*";

  procedure Collect(R1,R2: in Radius; DF: in Natural; F: in out RepMatrix; G: in out Components) is
--- collect higher order terms (not symmetric)
  begin
    if DF <= PDeg then return; end if;
    for D in reverse PDeg .. DF-1 loop
      UpMultAdd(R2,F(0,D+1),F(0,D));
      for I in 0 .. D loop
        UpMultAdd(R1,F(I+1,D-I),F(I,D-I));
      end loop;
    end loop;
    for I in Power loop
      Add(BallAt0(F(I,PDeg-I)),G(I,PDeg-I));
    end loop;
  end Collect;

  procedure MultAdd(DF,DG: in Natural; F,G: in Taylor2; H: in out Taylor2) is
--- only adds terms of order PDeg or less
    DF1:  constant Integer := Quasideg1(DF,F);
    DG1:  constant Integer := Quasideg1(DG,G);
    DFG:  constant Integer := IMin(DF+DG,PDeg);
    DFG1: constant Integer := IMin(DF1+DG1,PDeg);
    S: Scalar;
  begin
    for IH in reverse 0 .. DFG1 loop
      for JH in 0 .. DFG-IH loop
        S := ScalZero;
        for I in Imax(IH-DG1,0) .. IMin(IH,DF1) loop
          for J in Imax(JH-DG,0) .. IMin(JH,DF) loop
            MultAdd(F.C(I,J),G.C(IH-I,JH-J),S);
          end loop;
        end loop;
        Add(S,H.C(IH,JH));
      end loop;
    end loop;
  end MultAdd;

  procedure MultAdd(F,G: in Taylor2; H: in out Taylor2) is
    DF:   constant Integer := QuasiDeg(F);
    DG:   constant Integer := QuasiDeg(G);
    DFG:  constant Integer := DF+DG;
  begin
    if DF<0 or else DG<0 then return; end if;
    H.R := Min(Min(F.R,G.R),H.R);
    if Trunc or else DFG <= PDeg then          -- no higher order terms
      MultAdd(DF,DG,F,G,H);
      return;
    end if;
    declare
      DF1:  constant Integer := Quasideg1(DF,F);
      DG1:  constant Integer := Quasideg1(DG,G);
      DFG1: constant Integer := DF1+DG1;
      S: Scalar;
      Higher: RepMatrix := EZero;
    begin
      for IH in reverse 0 .. DFG1 loop
        for JH in 0 .. DFG-IH loop
          S := ScalZero;
          for I in Imax(IH-DG1,0) .. IMin(IH,DF1) loop
            for J in Imax(JH-DG,0) .. IMin(JH,DF) loop
              MultAdd(F.C(I,J),G.C(IH-I,JH-J),S);
            end loop;
          end loop;
          if IH+JH <= PDeg then
            Add(S,H.C(IH,JH));
          else
            Higher(IH,JH) := MaxNorm(S);
          end if;
        end loop;
      end loop;
      Collect(H.R(1),H.R(2),DFG,Higher,H.C);
    end;
  end MultAdd;

  procedure Scale(S1,S2: in Scalar; F: in out Taylor2) is
    D: constant Integer := QuasiDeg(F);
    T1,T2: Scalar;
  begin
    F.R(1) := Inf(Scal(F.R(1))/MaxNorm(S1));
    F.R(2) := Inf(Scal(F.R(2))/MaxNorm(S2));
    if D<1 then return; end if;
    T1 := ScalOne;
    for I in 0 .. D loop
      T2 := T1;
      for J in 0 .. D-I loop
        Mult(T2,F.C(I,J));
        Mult(S2,T2);
      end loop;
      Mult(S1,T1);
    end loop;
  end Scale;

  function Val(S1,S2: Scalar; F: Taylor2) return Scalar is
    R: constant RadPair :=  (MaxNorm(S1),MaxNorm(S2));
    S3,S4: Scalar;
  begin
    if not InDomain(R,F) then Error("Taylors2.Val domain error"); end if;
    S3 := ScalZero;
    for I in reverse Power loop
      S4 := ScalZero;
      for J in reverse 0 .. PDeg-I loop
        Mult(S2,S4);
        Add(F.C(I,J),S4);
      end loop;
      Mult(S1,S3);
      Add(S4,S3);
    end loop;
    return Disk(S3);
  end Val;

  procedure PolyDerU(F: in Taylor2; G: out Taylor2) is
  begin
    SetZero(F.R,G);
    for I in 1 .. PDeg loop
      for J in 0 .. PDeg-I loop
        Prod(I,F.C(I,J),G.C(I-1,J));
      end loop;
    end loop;
  end PolyDerU;

  procedure PolyDerV(F: in Taylor2; G: out Taylor2) is
  begin
    SetZero(F.R,G);
    for J in 1 .. PDeg loop
      for I in 0 .. PDeg-J loop
        Prod(J,F.C(I,J),G.C(I,J-1));
      end loop;
    end loop;
  end PolyDerV;

  function DerFac1(R1,R2: Radius) return Scalar is
--- norm factor for 1st derivative
--- R1: old domain, R2: new domain
    M: constant Scalar := Inv(Log(Scal(R1)/R2));
  begin
    if R1 <= R2 then
      Error("Taylors2.DerFac1 error");
    end if;
    if M<ScalOne then
      return ScalOne/R1;
    else
      return (M/R2)*Exp(-ScalOne);
    end if;
  end DerFac1;

  function DerFac2(R1,R2: Radius) return Scalar is
--- norm factor for 2nd derivative
--- R1: old domain, R2: new domain
    L: constant Scalar := Log(Scal(R1)/R2);
    M: constant Scalar := (L+ScalTwo+Sqrt(L*L+ScalFour))/(Two*L);
  begin
    if R1 <= R2 then
      Error("Taylors2.DerFac1 error");
    end if;
    if M<ScalTwo then
      return Two*Sqr(ScalOne/R1);
    else
      return (M/R2)*((M-ScalOne)/R2)*Exp(-(Two*M-ScalOne)/(M-ScalOne));
    end if;
  end DerFac2;

  procedure Ders(R: in RadPair; F0: in Taylor2; F: out Taylor2) is
  begin
    F := F0;
    SetRho(R,F);
  end Ders;

  procedure Ders(R: in RadPair; F0: in Taylor2; F,Fu,Fv: out Taylor2) is
  begin
    Ders(R,F0,F);
    PolyDerU(F,Fu);
    PolyDerV(F,Fv);
    if IsNumbers(F) then return; end if;
    declare
      U1: constant Scalar := DerFac1(F0.R(1),R(1));
      V1: constant Scalar := DerFac1(F0.R(2),R(2));
      G: Taylor2;
    begin
      ModNumbers(F,G);
      MultAdd(U1,G,Fu);
      MultAdd(V1,G,Fv);
    end;
  end Ders;

  procedure Ders(R: in RadPair; F0: in Taylor2; F,Fu,Fv,Fuu,Fuv,Fvv: out Taylor2) is
  begin
    Ders(R,F0,F);
    PolyDerU(F,Fu);
    PolyDerV(F,Fv);
    PolyDerU(Fu,Fuu);
    PolyDerV(Fu,Fuv);
    PolyDerV(Fv,Fvv);
    if IsNumbers(F) then return; end if;
    declare
      U1: constant Scalar := DerFac1(F0.R(1),R(1));
      V1: constant Scalar := DerFac1(F0.R(2),R(2));
      U2: constant Scalar := DerFac2(F0.R(1),R(1));
      V2: constant Scalar := DerFac2(F0.R(2),R(2));
      G: Taylor2;
    begin
      ModNumbers(Fu,G);
      MultAdd(V1,G,Fuv);
      MultAdd(Two*U1,G,Fuu);
      ModNumbers(Fv,G);
      MultAdd(U1,G,Fuv);
      MultAdd(Two*V1,G,Fvv);
      ModNumbers(F,G);
      MultAdd(U1,G,Fu);
      MultAdd(V1,G,Fv);
      MultAdd(U2,G,Fuu);
      MultAdd(V2,G,Fvv);
      MultAdd(U1*V1,G,Fuv);
    end;
  end Ders;

  procedure Ders(D: in Natural; R: in RadPair; F: in Taylor2; G: out Taylor2Mat) is
  begin
    if D=0 then
      Ders(R,F,G(0,0));
    elsif D=1 then
      Ders(R,F,G(0,0),G(1,0),G(0,1));
    elsif D=2 then
      Ders(R,F,G(0,0),G(1,0),G(0,1),G(2,0),G(1,1),G(0,2));
    else
      Error("Taylors2.Ders error: not implemented");
    end if;
  end Ders;

  procedure Show1(N: in String; H: in Taylor2) is
    W: constant Positive := 3;
    C: Scalar;
  begin
    if IsZero(H) then
      Put(N & "  Zero_Taylor2");
      New_Line;
    else
      Show2(N & " Rho",H.R(1),H.R(2));
      New_Line;
      for I in Power loop
        for J in 0 .. PDeg-I loop
          C := H.C(I,J);
          if C /= ScalZero then
            Show1(N & Strng(I,W) & Strng(J,W),C);
          end if;
        end loop;
      end loop;
    end if;
    New_Line;
  end Show1;

  procedure Show2(N: in String; H1,H2: in Taylor2) is
    W: constant Positive := 4;
    C1,C2: Scalar;
  begin
    if IsZero(H1) and IsZero(H2) then
      Put(N & "  Zero_Taylors2");
      New_Line;
    else
      Show2(N & " 1st Rho",H1.R(1),H2.R(1));
      Show2(N & " 2nd Rho",H1.R(2),H2.R(2));
      New_Line;
      for I in Power loop
        for J in 0 .. PDeg-I loop
          C1 := H1.C(I,J);
          C2 := H2.C(I,J);
          if not (IsZero(C1) and IsZero(C2)) then
            Show2(N & Strng(I,W) & Strng(J,W),C1,C2);
          end if;
        end loop;
      end loop;
    end if;
    New_Line;
  end Show2;

  procedure Get(F: in File_Type; H: out Taylor2; Decimal: in Boolean := False) is
    More: Boolean := False;
    D: Natural;
    C: Scalar;
  begin
    H := Tay_Zero;
    Get(F,H.R(1),Decimal);
    Get(F,H.R(2),Decimal);
    Get(F,D);
    D := D-1;
    declare
      Higher: RepMatrix := ZeroRepMat(D);
    begin
      for I in 0 .. D loop
        for J in 0 .. D-I loop
          Get(F,C,Decimal);
          if I+J <= PDeg then
            H.C(I,J) := C;
          else
            Higher(I,J) := MaxNorm(C);
            More := True;
          end if;
        end loop;
      end loop;
      if More and not Trunc then
        Collect(H.R(1),H.R(2),D,Higher,H.C);
      end if;
    end;
  end Get;

  procedure Put(F: in File_Type; H: in Taylor2; Decimal: in Boolean := False) is
    D: constant Integer := IMax(QuasiDeg(H),0);
  begin
    Put(F,H.R(1),Decimal);
    Put(F,H.R(2),Decimal);
    Put(F,D+1);
    for I in 0 .. D loop
      for J in 0 .. D-I loop
        Put(F,H.C(I,J),Decimal);
      end loop;
    end loop;
  end Put;

  procedure Read(Name: in String; H: out Taylor2; Decimal: in Boolean := False) is
    F: File_Type;
  begin
    Put_Line("Reading " & Name);
    Open(F,In_File,Name);
    Get(F,H,Decimal);
    Close(F);
  end Read;

  procedure Write(Name: in String; H: in Taylor2; Decimal: in Boolean := False) is
    F: File_Type;
  begin
    Put_Line("Writing " & Name);
    Create(F,Out_File,Name);
    Put(F,H,Decimal);
    Close(F);
  end Write;

begin

  Tay_Zero := Zero_Taylor2;

end Taylors2;