
package body Roots is

  use IR;

  procedure Decreasing(S1,S2: in out Scalar) is
  begin
    if MaxNorm(S1)<MaxNorm(S2) then Swap(S1,S2); end if;
  end Decreasing;

  procedure RRoots2(B,C: in Scalar; U1,U2,V: in out Scalar) is
  begin
    Quot(B,Flt(-2),U1);  -- U1 := -B/2   (fixed)
    Prod(U1,U1,U2);
    Sub(C,U2);           -- U2 := (B/2)**2-C
    if Sign(U2)<0 then
      Neg(U2);
      Sqrt(U2,V);        -- V := Sqrt(C-(B/2)**2)
      Copy(U1,U2);       -- U2 := -B/2
    else
      Sqrt(U2,V);        -- V := Sqrt((B/2)**2-C)
      Diff(U1,V,U2);     -- U2 := -B/2-V
      Add(V,U1);         -- U1 := -B/2+V
      Decreasing(U1,U2);
      Quot(C,U1,U2);     -- U2 := C/U1  (for precision)
      SetZero(V);
    end if;
  end RRoots2;

  procedure CRoots2(B,C: in Scalar; U1,U2,V: in out Scalar) is
  begin
    Quot(B,Flt(-2),U1);  -- U1 := -B/2   (fixed)
    Prod(U1,U1,U2);
    Sub(C,U2);           -- U2 := (B/2)**2-C
    Sqrt(U2,V);          -- V := Sqrt((B/2)**2-C)
    Diff(U1,V,U2);       -- U2 := -B/2-V
    Add(V,U1);           -- U1 := -B/2+V
    Decreasing(U1,U2);
    Quot(C,U1,U2);       -- U2 := C/U1  (for precision)
    SetZero(V);
  end CRoots2;

  procedure Roots2(B,C: in Scalar; U1,U2,V: in out Scalar) is
  begin
    if SComplex then
      CRoots2(B,C,U1,U2,V);
    else
      RRoots2(B,C,U1,U2,V);
    end if;
  end Roots2;

  procedure NewtonRoots2(B,C: in Scalar; U1,U2,V: in out Scalar; Steps: in Natural := 64) is
  begin
    if SComplex then
      CRoots2(B,C,U1,U2,V);
    else
      RRoots2(B,C,U1,U2,V);
    end if;
    if IsZero(V) then
      ImproveRoots2(B,C,U1,U2,Steps);
    end if;
  end NewtonRoots2;

  procedure Decreasing(S1,S2,S3: in out Scalar) is
  begin
    if MaxNorm(S2)<MaxNorm(S3) then Swap(S2,S3); end if;
    if MaxNorm(S1)<MaxNorm(S2) then Swap(S1,S2); end if;
    if MaxNorm(S2)<MaxNorm(S3) then Swap(S2,S3); end if;
  end Decreasing;

  procedure RRoots3(B,C,D: in Scalar; U0,U1,U2,V: in out Scalar) is
    Q:  constant Scalar := (B*B-Flt(3)*C)/Flt(9);
    R:  constant Scalar := ((Flt(2)*B*B-Flt(9)*C)*B+Flt(27)*D)/Flt(54);
    Q3: constant Scalar := Q*Q*Q;
    B3: constant Scalar := B/Three;
    AA,BB: Scalar;
  begin
    AA := R*R-Q3;
    if Sign(AA)<0 then
      AA := ArcCos(R/Sqrt(Q3));
      BB := Flt(4)*ArcCos(Scal(0));  -- 2*Pi
      V := Flt(-2)*Sqrt(Q);
      U0 := V*Cos(AA/Three)     -B3;
      U1 := V*Cos((AA+BB)/Three)-B3;
      U2 := V*Cos((AA-BB)/Three)-B3;
      Decreasing(U0,U1,U2);
      SetZero(V);
    else
      Copy(AA,BB);
      Sqrt(BB,AA);
      if Sign(R)<0 then
        Diff(AA,R,BB);
        Root(3,BB,AA);
      else
        Sum(AA,R,BB);
        Root(3,BB,AA);
        Neg(AA);
      end if;
      if IsZero(AA) then SetZero(BB); else Quot(Q,AA,BB); end if;
      Sum(AA,BB,V);
      Diff(V,B3,U0);
      Mult(NegHalf,V);
      Diff(V,B3,U1);
      Copy(U1,U2);
      V := Half*Sqrt(Scal(3))*(AA-BB);
    end if;
  end RRoots3;

  procedure CRoots3(B,C,D: in Scalar; U0,U1,U2,V: in out Scalar) is
    Q:  constant Scalar := (B*B-Flt(3)*C)/Flt(9);
    R:  constant Scalar := ((Flt(2)*B*B-Flt(9)*C)*B+Flt(27)*D)/Flt(54);
    Q3: constant Scalar := Q*Q*Q;
    B3: constant Scalar := B/Three;
    AA,BB: Scalar;
  begin
    AA := Sqrt(R*R-Q3);
    V := Adjoint(R)*AA;
    Real_Part(V);
    if Sign(V)<0 then Neg(AA); end if;
    Root(3,R+AA,AA);
    Neg(AA);
    if IsZero(AA) then SetZero(BB); else Quot(Q,AA,BB); end if;
    Sum(AA,BB,V);
    Diff(V,B3,U0);
    U1 := NegHalf*V-B3;
    V := Half*Sqrt(Scal(-3))*(AA-BB);
    Diff(U1,V,U2);
    Add(V,U1);
    Decreasing(U0,U1,U2);
    SetZero(V);
  end CRoots3;

  procedure Roots3(B,C,D: in Scalar; U0,U1,U2,V: in out Scalar) is
  begin
    if SComplex then
      CRoots3(B,C,D,U0,U1,U2,V);
    else
      RRoots3(B,C,D,U0,U1,U2,V);
    end if;
  end Roots3;

  procedure NewtonRoots3(B,C,D: in Scalar; U0,U1,U2,V: in out Scalar; Steps: in Natural := 64) is
  begin
    if SComplex then
      CRoots3(B,C,D,U0,U1,U2,V);
    else
      RRoots3(B,C,D,U0,U1,U2,V);
    end if;
    ImproveRoot3(B,C,D,U0,Steps);
    if IsZero(V) then
      ImproveRoots3(B,C,D,U1,U2,Steps);
    end if;
  end NewtonRoots3;

end Roots;
