end;

    { Double QR step involving rows L to En and columns M to En }
    for K := M to Na do
      begin
        NotLas := (K <> Na);
        if (K = M) then goto 170;
        P := H[K,K - 1];
        Q := H[K + 1,K - 1];
        R := 0.0;
        if NotLas then R := H[K + 2,K - 1];
        X := Abs(P) + Abs(Q) + Abs(R);
        if X = 0.0 then goto 260;
        P := P / X;
        Q := Q / X;
        R := R / X;
170:    S := DSgn(Sqrt(P * P + Q * Q + R * R), P);
        if K <> M then
          H[K,K - 1] := - S * X
        else if L <> M then
          H[K,K - 1] := - H[K,K - 1];
        P := P + S;
        X := P / S;
        Y := Q / S;
        Zz := R / S;
        Q := Q / P;
        R := R / P;
        if NotLas then goto 225;

        { Row modification }
        for J := K to Ubound do
          begin
            P := H[K,J] + Q * H[K + 1,J];
            H[K,J] := H[K,J] - P * X;
            H[K + 1,J] := H[K + 1,J] - P * Y;
          end;

        J := Min(En, K + 3);

        { Column modification }
        for I := Lbound to J do
          begin
            P := X * H[I,K] + Y * H[I,K + 1];
            H[I,K] := H[I,K] - P;
            H[I,K + 1] := H[I,K + 1] - P * Q;
          end;

        { Accumulate transformations }
        for I := I_low to I_igh do
          begin
            P := X * Z[I,K] + Y * Z[I,K + 1];
            Z[I,K] := Z[I,K] - P;
            Z[I,K + 1] := Z[I,K + 1] - P * Q;
          end;
        goto 260;

225:
        { Row modification }
        for J := K to Ubound do
          begin
            P := H[K,J] + Q * H[K + 1,J] + R * H[K + 2,J];
            H[K,J] := H[K,J] - P * X;
            H[K + 1,J] := H[K + 1,J] - P * Y;
            H[K + 2,J] := H[K + 2,J] - P * Zz;
          end;

        J := Min(En, K + 3);

        { Column modification }
        for I := Lbound to J do
          begin
            P := X * H[I,K] + Y * H[I,K + 1] + Zz * H[I,K + 2];
            H[I,K] := H[I,K] - P;
            H[I,K + 1] := H[I,K + 1] - P * Q;
            H[I,K + 2] := H[I,K + 2] - P * R;
          end;

        { Accumulate transformations }
        for I := I_low to I_igh do
          begin
            P := X * Z[I,K] + Y * Z[I,K + 1] + Zz * Z[I,K + 2];
            Z[I,K] := Z[I,K] - P;
            Z[I,K + 1] := Z[I,K + 1] - P * Q;
            Z[I,K + 2] := Z[I,K + 2] - P * R;
          end;

260:  end;

    goto 70;

270: { One root found }
    H[En,En] := X + T;
    Wr[En] := H[En,En];
    Wi[En] := 0.0;
    En := Na;
    goto 60;

280: { Two roots found }
    P := 0.5 * (Y - X);
    Q := P * P + W;
    Zz := Sqrt(Abs(Q));
    H[En,En] := X + T;
    X := H[En,En];
    H[Na,Na] := Y + T;
    if Q < 0.0 then goto 320;

    { Real pair }
    Zz := P + DSgn(Zz, P);
    Wr[Na] := X + Zz;
    Wr[En] := Wr[Na];
    if Zz <> 0.0 then Wr[En] := X - W / Zz;
    Wi[Na] := 0.0;
    Wi[En] := 0.0;
    X := H[En,Na];
    S := Abs(X) + Abs(Zz);
    P := X / S;
    Q := Zz / S;
    R := Sqrt(P * P + Q * Q);
    P := P / R;
    Q := Q / R;

    { Row modification }
    for J := Na to Ubound do
      begin
        Zz := H[Na,J];
        H[Na,J] := Q * Zz + P * H[En,J];
        H[En,J] := Q * H[En,J] - P * Zz;
      end;

    { Column modification }
    for I := Lbound to En do
      begin
        Zz := H[I,Na];
        H[I,Na] := Q * Zz + P * H[I,En];
        H[I,En] := Q * H[I,En] - P * Zz;
      end;

    { Accumulate transformations }
    for I := I_low to I_igh do
      begin
        Zz := Z[I,Na];
        Z[I,Na] := Q * Zz + P * Z[I,En];
        Z[I,En] := Q * Z[I,En] - P * Zz;
      end;

    goto 330;

320: { Complex pair }
    Wr[Na] := X + P;
    Wr[En] := Wr[Na];
    Wi[Na] := Zz;
    Wi[En] := - Zz;

330:
    En := Enm2;
    goto 60;

340:
    if Norm = 0.0 then Exit;

    { All roots found. Backsubstitute to find
      vectors of upper triangular form }
    for En := Ubound downto Lbound do
      begin
        P := Wr[En];
        Q := Wi[En];
        Na := En - 1;
        if Q < 0.0 then
          goto 710
        else if Q = 0.0 then
          goto 600
        else
          goto 800;

600:    { Real vector }
        M := En;
        H[En,En] := 1.0;
        if Na < Lbound then goto 800;

        for I := Na downto Lbound do
          begin
            W := H[I,I] - P;
            R := 0.0;

            for J := M to En do
              R := R + H[I,J] * H[J,En];

            if Wi[I] >= 0.0 then goto 630;
            Zz := W;
            S := R;
            goto 700;
630:        M := I;
            if Wi[I] <> 0.0 then goto 640;
            T := W;
            if T <> 0.0 then goto 635;
            Tst1 := Norm;
            T := Tst1;
            repeat
              T := 0.01 * T;
              Tst2 := Norm + T;
            until Tst2 <= Tst1;
635:        H[I,En] := - R / T;
            goto 680;

640:        { Solve real equations }
            X := H[I,I + 1];
            Y := H[I + 1,I];
            Q := (Wr[I] - P) * (Wr[I] - P) + Wi[I] * Wi[I];
            T := (X * S - Zz * R) / Q;
            H[I,En] := T;
            if Abs(X) > Abs(Zz) then
              H[I + 1,En] := (- R - W * T) / X
            else
              H[I + 1,En] := (- S - Y * T) / Zz;

680:        { Overflow control }
            T := Abs(H[I,En]);
            if T = 0.0 then goto 700;
            Tst1 := T;
            Tst2 := Tst1 + 1.0 / Tst1;
            if Tst2 > Tst1 then goto 700;
            for J := I to En do
              H[J,En] := H[J,En] / T;
700:      end;
        { End real vector }
        goto 800;

        { Complex vector }
710:    M := Na;

        { Last vector component chosen imaginary so that
          eigenvector matrix is triangular }
        if Abs(H[En,Na]) > Abs(H[Na,En]) then
          begin
            H[Na,Na] := Q / H[En,Na];
            H[Na,En] := - (H[En,En] - P) / H[En,Na];
          end
        else
          Cdiv(0.0, - H[Na,En], H[Na,Na] - P, Q, H[Na,Na], H[Na,En]);

        H[En,Na] := 0.0;
        H[En,En] := 1.0;
        Enm2 := Na - 1;
        if Enm2 < Lbound then goto 800;

        for I := Enm2 downto Lbound do
          begin
            W := H[I,I] - P;
            Ra := 0.0;
            Sa := 0.0;

            for J := M to En do
              begin
                Ra := Ra + H[I,J] * H[J,Na];
                Sa := Sa + H[I,J] * H[J,En];
              end;

            if Wi[I] >= 0.0 then goto 770;
            Zz := W;
            R := Ra;
            S := Sa;
            goto 795;
770:        M := I;
            if Wi[I] <> 0.0 then goto 780;
            Cdiv(- Ra, - Sa, W, Q, H[I,Na], H[I,En]);
            goto 790;

            { Solve complex equations }
780:        X := H[I,I + 1];
            Y := H[I + 1,I];
            Vr := (Wr[I] - P) * (Wr[I] - P) + Wi[I] * Wi[I] - Q * Q;
            Vi := (Wr[I] - P) * 2.0 * Q;
            if (Vr = 0.0) and (Vi = 0.0) then
              begin
                Tst1 := Norm * (Abs(W) + Abs(Q) + Abs(X) + Abs(Y) + Abs(Zz));
                Vr := Tst1;
                repeat
                  Vr := 0.01 * Vr;
                  Tst2 := Tst1 + Vr;
                until Tst2 <= Tst1;
              end;
            Cdiv(X * R - Zz * Ra + Q * Sa, X * S - Zz * Sa - Q * Ra, Vr, Vi, H[I,Na], H[I,En]);
            if Abs(X) > Abs(Zz) + Abs(Q) then
              begin
                H[I + 1,Na] := (- Ra - W * H[I,Na] + Q * H[I,En]) / X;
                H[I + 1,En] := (- Sa - W * H[I,En] - Q * H[I,Na]) / X;
              end
            else
              Cdiv(- R - Y * H[I,Na], - S - Y * H[I,En], Zz, Q, H[I + 1,Na], H[I + 1,En]);

790:        { Overflow control }
            T := Max(Abs(H[I,Na]), Abs(H[I,En]));
            if T = 0.0 then goto 795;
            Tst1 := T;
            Tst2 := Tst1 + 1.0 / Tst1;
            if Tst2 > Tst1 then goto 795;
            for J := I to En do
              begin
                H[J,Na] := H[J,Na] / T;
                H[J,En] := H[J,En] / T;
              end;

795:      end;
      { End complex vector }
800:  end;

    { End back substitution.
      Vectors of isolated roots }
    for I := Lbound to Ubound do
      if (I < I_low) or (I > I_igh) then
        for J := I to Ubound do
          Z[I,J] := H[I,J];

    { Multiply by transformation matrix to give
      vectors of original full matrix. }
    for J := Ubound downto I_low do
      begin
        M := Min(J, I_igh);
        for I := I_low to I_igh do
          begin
            Zz := 0.0;
            for K := I_low to M do
              Zz := Zz + Z[I,K] * H[K,J];
            Z[I,J] := Zz;
          end;
      end;

    Hqr2 := 0;
  end;

  procedure BalBak(Z : TMatrix; Lbound, Ubound, I_low, I_igh : Integer;
                   Scale : TVector; M : Integer);
{ ----------------------------------------------------------------------
  This procedure is a translation of the EISPACK subroutine Balbak

  This procedure forms the eigenvectors of a real general matrix
  by back transforming those of the corresponding balanced matrix
  determined by Balance.

  On input:

    Z contains the real and imaginary parts of the eigenvectors
    to be back transformed.

    Lbound, Ubound are the lowest and highest indices
    of the elements of Z

    I_low and I_igh are integers determined by Balance.

    Scale contains information determining the permutations
    and scaling factors used by Balance.

    M is the index of the latest column of Z to be back transformed.

  On output:

    Z contains the real and imaginary parts of the transformed
    eigenvectors in its columns Lbound..M
  ---------------------------------------------------------------------- }
  var
    I, J, K : Integer;
    S : Float;
  begin
    if M < Lbound then Exit;

    if I_igh <> I_low then
      for I := I_low to I_igh do
        begin
          S := Scale[I];
          { Left hand eigenvectors are back transformed if the
            foregoing statement is replaced by S := 1.0 / Scale[I] }
          for J := Lbound to M do
            Z[I,J] := Z[I,J] * S;
        end;

    for I := (I_low - 1) downto Lbound do
      begin
        K := Round(Scale[I]);
        if K <> I then
          for J := Lbound to M do
            Swap(Z[I,J], Z[K,J]);
      end;

    for I := (I_igh + 1) to Ubound do
      begin
        K := Round(Scale[I]);
        if K <> I then
          for J := Lbound to M do
            Swap(Z[I,J], Z[K,J]);
      end;
  end;

  function EigenVals(A : TMatrix; Lbound, Ubound : Integer;
                     Lambda_Re, Lambda_Im : TVector) : Integer;
  var
    I_low, I_igh : Integer;
    Scale : TVector;
    I_int : TIntVector;
  begin
    DimVector(Scale, Ubound);
    DimVector(I_Int, Ubound);

    Balance(A, Lbound, Ubound, I_low, I_igh, Scale);
    ElmHes(A, Lbound, Ubound, I_low, I_igh, I_int);
    EigenVals := Hqr(A, Lbound, Ubound, I_low, I_igh, Lambda_Re, Lambda_Im);
  end;

  function EigenVect(A : TMatrix; Lbound, Ubound : Integer;
                     Lambda_Re, Lambda_Im : TVector; V : TMatrix) : Integer;
  var
    I_low, I_igh, ErrCode : Integer;
    Scale : TVector;
    I_Int : TIntVector;
  begin
    DimVector(Scale, Ubound);
    DimVector(I_Int, Ubound);

    Balance(A, Lbound, Ubound, I_low, I_igh, Scale);
    ElmHes(A, Lbound, Ubound, I_low, I_igh, I_int);
    Eltran(A, Lbound, Ubound, I_low, I_igh, I_int, V);
    ErrCode := Hqr2(A, Lbound, Ubound, I_low, I_igh, Lambda_Re, Lambda_Im, V);
    if ErrCode = 0 then BalBak(V, Lbound, Ubound, I_low, I_igh, Scale, Ubound);

    EigenVect := ErrCode;
  end;

  procedure DivLargest(V : TVector; Lbound, Ubound : Integer;
                       var Largest : Float);
  var
    I : Integer;
  begin
    Largest := V[Lbound];
    for I := Succ(Lbound) to Ubound do
      if Abs(V[I]) > Abs(Largest) then
        Largest := V[I];
    for I := Lbound to Ubound do
      V[I] := V[I] / Largest;
  end;

  function RootPol(Coef : TVector; Deg : Integer;
                   Xr, Xi : TVector) : Integer;
  var
    A            : TMatrix;  { Companion matrix }
    N            : Integer;  { Size of matrix }
    I_low, I_igh : Integer;  { Used by Balance }
    Scale        : TVector;  { Used by Balance }
    Nr           : Integer;  { Number of real roots }
    I, J         : Integer;  { Loop variables }
    ErrCode      : Integer;  { Error code }
  begin
    N := Pred(Deg);
    DimMatrix(A, N, N);
    DimVector(Scale, N);

    { Set up the companion matrix (to save space, begin at index 0) }
    for J := 0 to N do
      A[0,J] := - Coef[Deg - J - 1] / Coef[Deg];
    for J := 0 to Pred(N) do
      A[J + 1,J] := 1.0;

    { The roots of the polynomial are the eigenvalues of the companion matrix }
    Balance(A, 0, N, I_low, I_igh, Scale);
    ErrCode := Hqr(A, 0, N, I_low, I_igh, Xr, Xi);

    if ErrCode <> 0 then
      begin
        RootPol := ErrCode;
        Exit;
      end;

    { Count real roots }
    Nr := 0;
    for I := 0 to N do
      if Xi[I] = 0.0 then
        Inc(Nr);

    { Transfer roots from 0..(Deg - 1) to 1..Deg }
    for I := N downto 0 do
      begin
        J := I + 1;
        Xr[J] := Xr[I];
        Xi[J] := Xi[I];
      end;

    RootPol := Nr;
  end;

end.