File : funs2-num.adb

with Messages;
use Messages;

package body Funs2.Num is

  function NumInv(F: Fun) return Fun is
    Iter: constant Integer := 15;
    Eps: constant Radius := 5.0E-16;
    Err: Rep;
    G,H: Fun;
  begin
    G.A := F.A;
    G.E := F.E;
    if IsZero(F.Q) then
      G.P := NumInv(F.P);
      G.Q := TayZero;
    else
      Prod(F.P,F.P,G.P);
      Prod(F.Q,F.Q,G.Q);
      Mult(CoordT2(F.A),G.Q);
      G.Q := NumInv(G.P-G.Q);
      G.P := F.P*G.Q;
      G.Q := -F.Q*G.Q;
    end if;
    for I in 1 .. Iter loop    -- now iterate G -> (2-FG)G
      Prod(F,G,H);
      AddConst(-ScalOne,H);
      Err := Sup(Norm(H));
      AddConst(-ScalOne,H);
      Mult(H,G);
      Neg(G);
      if Err<Eps then return G; end if;
    end loop;
    Show("Funs2.Num.NumInv: precision only",Err,1);
    return G;
  end NumInv;

  procedure NumInv(L: in out Flop) is
    H: Fun;
  begin
    H := L(1)*Reflect(L(1))-L(2)*Reflect(L(2));
    H := NumInv(H);
    L(1) := H*Reflect(L(1));
    L(2) := -H*L(2);
  end NumInv;

  procedure VApprox(D: in Domain; V: out Fun) is
  begin
    if not VHave then
      Read("vapprox.data",VSave);
      VHave := True;
    end if;
    if VSave.A /= D.A then
      if Verbosity>2 then Show2("converting ",VSave.A,D.A); end if;
      ChangeVar(D.A,VSave);
    end if;
    V := VSave;
    SetRho(D.R,V);
  end VApprox;

  procedure NumCompZero(F,G: in Fun; R: in RadPair; V: in out Fun) is
--- solves G(x,V)+F(V,W)=0 by Newton's method
--- V can be zero but needs valid V.A
--- using also VApprox if initial V is zero or inaccurate
    OldVerbosity: constant Natural := Verbosity;
    TmpVerbosity: constant Natural := IMin(2,Verbosity);
    Iter: constant Integer := 30;
    N:    constant Integer := 4;
    Eps1: constant Radius  := 1.0E-05;
    Eps2: constant Radius  := 1.0E-08;
    Eps3: constant Radius  := 1.0E-15;
    DV:   constant Domain  := (V.A,R);
    Skip: Boolean := False;
    Err: Rep;
    FX,FY,GX,GY,FGV: Fun;
    L: Flop;
  begin
    TraceEnter("Funs.Num.NumCompZero");
    Show("Computing approx V",3);
    Verbosity := TmpVerbosity;           --- too many domain violations
    SetRho(R,V);
    if not IsZero(V) then
      FGV := CompY(V,G)+CompVW(V,F);
      Err := Sup(Norm(FGV));
      Skip := Err<Eps1;
    end if;
    if not Skip then VApprox(DV,V); end if;
    for I in 1 .. Iter loop
      if not Skip then
        FGV := CompY(V,G)+CompVW(V,F);
        Err := Sup(Norm(FGV));
      end if;
      Skip := False;
      if Err>Eps2 or else (I mod N)=1 then
        Der(R,F,FX,FY);
        Der(R,G,GX,GY);
        L(1) := CompY(V,GY)+CompVW(V,FX);
        L(2) := -CompVW(V,FY);
        NumInv(L);
      end if;
      V := V-L*FGV;
      exit when Err<Eps3;
    end loop;
    if Err>Eps3 then
      Show("Funs2.Num.NumCompZero: precision only",Err,1);
    end if;
    Verbosity := OldVerbosity;
    TraceLeave;
  end NumCompZero;

  procedure NumCompZeroD(F,G: in Fun; R: RadPair; V: in out Fun) is
--- NumCompZero for derivatives
    FX,GY,H: Fun;
  begin
    Der(Rho(F),F,FX,H);
    Der(Rho(G),G,H,GY);
    NumCompZero(FX,GY,R,V);
  end NumCompZeroD;

  function GuessRho(F: Fun) return RadPair is
    D0: constant Integer := PDeg/5;
    D1: constant Integer := PDeg-D0;
  begin
    return Min(GuessRho(D0,D1,F.P),GuessRho(D0,D1,F.Q));
  end GuessRho;

  procedure RhoInfo(N: in String; F: in Fun; Rn: in RadPair) is
    Rh: constant RadPair := Rho(F);
    Rg: constant RadPair := GuessRho(F);
  begin
    Show1(N,F.A);
    Show(N & " guess radii ",Rg(1),Rg(2));
    Show(N & "  have radii ",Rh(1),Rh(2));
    Show(N & "  need radii ",Rn(1),Rn(2));
    Show(N & " need / have ",Rn(1)/Rh(1),Rn(2)/Rh(2));
    New_Line;
  end RhoInfo;

end Funs2.Num;