with Strings;
use Strings;

package body Bisection is

  procedure FindSigns(X0,X: in out Scalar; S0,S: out Integer) is
  begin
    S0 := SignF(X0);
    if S0=0 then
      Copy(X0,X);
      S := S0;
      return;
    end if;
    S := SignF(X);
    if S /= S0 then return; end if;
    declare
      KMax: constant Positive := 2**NSign;
      DK: Positive := KMax;
      K: Positive;
      DX: Scalar;
    begin
      Diff(X,X0,DX);
      Mult(One/Flt(KMax),DX);     -- DX := (X-X0)/Flt(KMax)
      while DK>1 loop
        K := DK/2;
        while K<KMax loop
          Copy(DX,X);
          Mult(Flt(K),X);
          Add(X0,X);              -- X := X0+K*DX
          S := SignF(X);
          if S /= S0 then
            Mult(Flt(K-DK/2),DX);
            Add(DX,X0);           -- X0 := X0+(K-DK/2)*DX
            return;
          end if;
          K := K+DK;
        end loop;
        DK := DK/2;
      end loop;
    end; --- no sign change found
  end FindSigns;

  procedure Find_Signs(X1,X2: in out Flt; Y1,Y2: out Flt) is
    Done: exception;
    pragma Warnings (Off);
  begin
    Y1 := F(X1);
    Y2 := F(X2);
    if Y1*Y2 <= Zero then return; end if;
    declare
      Eps: constant Flt := Half;
      Neg: constant Boolean := (Y1<Zero);
      X,Xt,Y,Yt: Flt;

      function G(S: Flt) return Flt is
      begin
        if Neg then return -F(S); else return F(S); end if;
      end G;

      procedure Find_Smaller is
        ---     expecting 0<G(X1)<G(X2) and G'(X1)*DX<0
        DX: Flt := Eps*(X2-X1);
      begin
        X := X1;
        loop
          Xt := X;
          X := X1+DX;
          exit when X=Xt;
          Y := G(X);
          if Y<Y1 then return; end if;
          DX := Half*DX;
        end loop;
        Show0("Bisection.Find_Signs: no smaller G-value found");
        raise Done;
      end Find_Smaller;

      procedure Find_Neg is
        ----    expecting G(X)<G(X1)<G(X2)
        Fac1: constant Flt := Flt(5063)/Flt(8192);
        Fac2: constant Flt := Flt(3129)/Flt(8192);
      begin
        loop
          if abs(X-X1)<abs(X-X2) then -- X closer to X1
            Xt := Fac2*X+Fac1*X2;
            exit when Xt=X;
            Yt := G(Xt);
            if Yt<Y then
              if Yt<Zero then
                X1 := X; Y1 := Y;
                X2 := Xt; Y2 := Yt;
                raise Done;          --- creates warning
              end if;
              X1 := X; Y1 := Y;
              X := Xt; Y := Yt;
            else
              X2 := Xt; Y2 := Yt;
            end if;
          else                        -- X closer to X2
            Xt := Fac2*X+Fac1*X1;
            exit when Xt=X;
            Yt := G(Xt);
            if Yt<Y then
              if Yt<Zero then
                X1 := X; Y1 := Y;
                X2 := Xt; Y2 := Yt;
                raise Done;
              end if;
              X2 := X; Y2 := Y;
              X := Xt; Y := Yt;
            else
              X1 := Xt; Y1 := Yt;
            end if;
          end if;
        end loop;
        Show0("Bisection.Find_Signs: no negative G-value found");
        raise Done;
      end Find_Neg;

    begin
      if Neg then Y1 := -Y1; Y2 := -Y2; end if;
      if Y2<Y1 then
        Xt := X2; X2 := X1; X1 := Xt;
        Yt := Y2; Y2 := Y1; Y1 := Yt;
      end if;
      Find_Smaller; --- now G(X)<G(X1)<G(X2)
      if Y<Zero then
        X2 := X; Y2 := Y;
        raise Done;
      end if;
      Find_Neg;
    exception
      when Done => if Neg then Y1 := -Y1; Y2 := -Y2; end if;
    end;
  end Find_Signs;

  procedure FindZero(XN,XP: in out Scalar; Err: in out Radius) is
    --- notice that (XP+XN)/2 could become XP or XN
    MaxCount: constant Natural := 4; -- max number of consecutive equal F-values
    Eps: constant Radius := Err;
    Count: Natural := 0;
    Y: Flt;
    YSave: Flt := Zero;
    X: Scalar;
  begin
    for N in 1 .. NFind loop
      Sum(XP,XN,X);
      Mult(Half,X);  -- X := (XP+XN)/2
      Y := F(X);
      exit when abs(Y)<=Eps;
      if Y=YSave then
        Count := Count+1;
        exit when Count=MaxCount;
      else
        Count := 0;
        YSave := Y;
      end if;
      if Y>Zero then Copy(X,XP); else Copy(X,XN); end if;
    end loop;
    Err := abs(Y);
  end FindZero;

end Bisection;
