File : zeros.adb

with Messages;
use Messages;

package body Zeros is

  Iter: constant Integer := 30;

  function FindZero(S0,S1,S2: Number; Eps1,Eps2: Radius) return Number is
    FoundPos,FoundNeg: Boolean := False;
    N: Integer := 8;
    S: Number := S0;
    Err: Rep;
    SP,SN,DS,DFS,Rad,K: Number;
  begin

    if not Z_Trunc then                     --- just add error estimate to the guess S0
      S := Center(S);
      DFS := DF(S);
      DS := F(S)/DFS;
      Rad := Rep(N+1)*Norm(DS)/Rep(N);      --- radius of ball
      S := S+DiskAt0(Rad);                  --- ball
      K := Scal_One-DF(S)/DFS;              --- contraction factor
      if MaxNorm(Rep(N+1)*K)>One then
        Error("Zeros.FindZero error: no contraction");
      end if;
      return S;
     end if;

    if S=Scal_Zero then                      --- bisection
      for I in 1 .. Iter loop
        if (I mod 2)=0 then N := I/2; else N := -(I/2); end if;
        S := Scal(1-N)*S1+Scal(N)*S2;
        if F(S) >= Scal_Zero then
          FoundPos := True;
          SP := S;
          exit when FoundNeg;
        else
          FoundNeg := True;
          SN := S;
          exit when FoundPos;
        end if;
      end loop;
      if not (FoundPos and FoundNeg) then
        Error("Zeros.FindZero error: no zero found");
      end if;
      for I in 1 .. Iter loop
        Err := MaxNorm(SP-SN);
        S := (SP+SN)/Scal_Two;
        exit when Err<Eps1;
        if F(S) >= Scal_Zero then SP := S; else SN := S; end if;
        if F(S) >= Scal_Zero then SP := S; else SN := S; end if;
      end loop;
      if Err>Eps1 then
        Show("Zeros.FindZero: bisection precision only",Err,1);
      end if;
    end if;

    for I in 1 .. Iter loop                 --- Newton
      DS := F(S)/DF(S);
      Err := MaxNorm(DS);
      S := S-DS;
      exit when Err<Eps2;
    end loop;
    if Err>Eps2 then
      Show("Zeros.FindZero: Newton precision only",Err,1);
    end if;
    return S;

  end FindZero;

end Zeros;