with Strings;
use Strings;

pragma Elaborate_all (Strings);

package body Newton is

  procedure NumFind0(S: in out Scalar; Steps: in Natural := 64; Verbosity: in Natural := 0) is
    I: Integer := 0;
    Eps0,Eps1: Flt := Flt'Last;
    FS,DS : Scalar;
    DFS: Operator;

    procedure ShowEps is
    begin
      if Steps>0 then
        Show1("FindZero " & Strng(I) & " Eps =",Eps1);
        if I=Steps then
          Show0("FindZero: error still decreasing");
        end if;
      end if;
    end ShowEps;

  begin
    if Numeric then
      while I<Steps loop       -- S -> S-F(S)/DF(S);
        I := I+1;
        Eps0 := Eps1;
        TF(S,FS,DFS);
        Quot(FS,DFS,DS);
        Eps1 := MaxNorm(DS);
        if Eps1 <= MinErr then
          exit when Eps1=Zero or else (Eps1 >= Eps0);
        end if;
        if Verbosity>4 then ShowEps; end if;
        Sub(DS,S);
      end loop;
    else
      Center(S);
      while I<Steps loop          -- S -> S-F(S)/DF(S);
        I := I+1;
        Eps0 := Eps1;
        TF(S,FS,DFS);
        Center(FS);
        Center(DFS);
        Quot(FS,DFS,DS);
        Eps1 := MaxNorm(DS);
        if Eps1 <= MinErr then
          exit when Eps1=Zero or else (Eps1 >= Eps0);
        end if;
        if Verbosity>4 then ShowEps; end if;
        Sub(DS,S);
        Center(S);
      end loop;
    end if;
    if Verbosity>3 then ShowEps; end if;
  end NumFind0;

  ------------------------------------------------

  procedure NumFindZero(S: in out Scalar; Steps: in Natural := 64; Verbosity: in Natural := 0) is

    procedure TF(X: in Scalar; Y: in out Scalar; L: in out Operator) is
    begin
      F(X,Y);
      DF(X,L);
    end TF;

    procedure Find is new NumFind0 (TF => TF);
  begin
    Find(S,Steps,Verbosity);
  end NumFindZero;

  procedure Find0(S: in out Scalar; Steps: in Natural := 64; Verbosity: in Natural := 0) is
  begin
    if Steps>0 then
      declare
        procedure Find is new NumFind0 (TF => TF);
      begin
        Find(S,Steps,Verbosity);
      end;
    end if;
    if not Numeric then
      declare
        N: constant Positive := 8; -- contraction factor is 1/(N+1)
        R,K: Flt;
        FS,DS,E: Scalar;
        DFS,L,Tmp: Operator;
      begin
        Center(S);
        TF(S,FS,DFS);
        Quot(FS,DFS,DS);               -- DS := F(S)/DF(S)
        R := MaxNorm(DS)*Flt(N+1)/Flt(N);
        BallAt0(R,E);                  -- Ball of radius ||DS||*(N+1)/N
        Add(E,S);                      -- S+Ball
        TF(S,FS,Tmp);
        Quot(Tmp,DFS,L);
        Add(-1,L);
        K := MaxNorm(L)*Flt(N+1);      -- K := ||DF(S+Ball)/DF(S)-1||*(N+1)
        if Verbosity>6 then
          Show1("K =",K);
        end if;
        if K>One then
          Show1("K =",K);
          raise No_Contraction;
        end if;
      end;
    end if;
  end Find0;

  procedure FindZero(S: in out Scalar; Steps: in Natural := 64; Verbosity: in Natural := 0) is

    procedure TF(X: in Scalar; Y: in out Scalar; L: in out Operator) is
    begin
      F(X,Y);
      DF(X,L);
    end TF;

    procedure Find is new Find0 (TF => TF);
  begin
    Find(S,Steps,Verbosity);
  end FindZero;

  procedure Improve0(S: in out Scalar; Steps: in Natural := 64) is
    --- S is already known to enclose the zero
    Eps0,Eps1: Flt := Flt'Last;
    FS,DS : Scalar;
    DFS: Operator;
  begin
    for I in 1 .. Steps loop  -- S -> S-F(S)/DF(S);
      Eps0 := Eps1;
      TF(S,FS,DFS);
      Quot(FS,DFS,DS);
      Eps1 := MaxNorm(DS);
      exit when Eps1=Zero or else (Eps1 >= Eps0);
      Sub(DS,S);
    end loop;
  end Improve0;

  procedure ImproveZero(S: in out Scalar; Steps: in Natural := 64) is

    procedure TF(X: in Scalar; Y: in out Scalar; L: in out Operator) is
    begin
      F(X,Y);
      DF(X,L);
    end TF;

    procedure Improve is new Improve0 (TF => TF);
  begin
    Improve(S,Steps);
  end ImproveZero;

  procedure InvNewton(L1: in Operator; L2: in out Operator; Steps: in Natural := 64; Verbosity: in Natural := 0) is
    N: constant Positive := 8; -- contraction factor is 1/(N+1)
    I: Integer := 0;
    Eps0,Eps1: Flt := Flt'Last;
    Del,Eps,R,K: Flt;
    L3,Tmp: Operator;

    procedure ShowEps is
    begin
      if Steps>0 then
        Show1("InvNewton" & I'Img & " Eps =",Eps1);
        if I=Steps then
          Show0("InvNewton: error still decreasing");
        end if;
      end if;
    end ShowEps;

  begin
    if Numeric then
      while I<Steps loop          -- L2 -> L2*(2-L1*L2)
        I := I+1;
        Eps0 := Eps1;
        Copy(L2,L3);
        Prod(L1,L2,Tmp);
        Neg(Tmp);
        Add(1,Tmp);
        Eps1 := MaxNorm(Tmp);
        exit when Eps1=Zero or else (Eps1 >= Eps0);
        if Verbosity>4 then ShowEps; end if;
        Add(1,Tmp);
        Prod(L3,Tmp,L2);
      end loop;
      if Verbosity>3 then ShowEps; end if;
    else
      Center(L2);
      while I<Steps loop          -- L2 -> L2*(2-L1*L2)
        I := I+1;
        Eps0 := Eps1;
        Copy(L2,L3);
        Prod(L1,L2,Tmp);
        Center(Tmp);
        Neg(Tmp);
        Add(1,Tmp);
        exit when Eps1=Zero or else (Eps1 >= Eps0);
        if Verbosity>4 then ShowEps; end if;
        Add(1,Tmp);
        Prod(L3,Tmp,L2);
        Center(L2);
      end loop;
      if Verbosity>3 then ShowEps; end if;
      Prod(L1,L2,L3);
      Add(-1,L3);
      Del := MaxNorm(L3);              -- Del := ||L1*L2-1||
      Eps := Del*MaxNorm(L2);          -- Eps := Del*||L2||
      R := Eps*Flt(N+1)/Flt(N);        -- R := Eps*(N+1)/N
      K := Two*(Del+R*MaxNorm(L1));    -- K := 2*Del+2*R*||P1||
      K := Flt(N+1)*K;
      if Verbosity>6 then
        Show1("K =",K);
      end if;
      if K>One then
        Show1("K =",K);
        raise No_Contraction;
      end if;
      BallAt0(R,Tmp);
      Add(Tmp,L2);
    end if;
  end InvNewton;

end Newton;
