"The Useful Program: define p0 and let it go (just make sure p0 has \
a positive root)";
"Note: when Mathematica finds roots of a polynomial, it puts them \
in this order: first the real roots, in increasing order, then the imaginary \
roots. This is why the adj subroutine works to get Mathematica to look for \
exactly the greatest positive real root.";

p0 = x^3 - 5  x - 1;

Plot[p0, {x, -5, 5}];

k = p0;
n = 10;
d = Exponent[p0, x];
adj = 
        Catch[Do[
            If[Im[Root[p0, d - i]] == 0, Throw[i]], {i, 0, 
              d - 1}]];
a0 = Floor[Root[p0, d - adj]];
Coeffs = 
      Table[Floor[
          Root[ k = 
              Simplify[(-x^d) (k  /.  
                      x -> (1/x + Floor[Root[k, d - adj]]))] , 
            d - adj]], {n}]
Frac = 
      a0 + (FromContinuedFraction[
              Coeffs])^(-1)
N[Frac]
N[Root[p0, d - adj]]

"Now, some steps, examples, and variations...";

"The original program";

p0 = k = x^2 - 5  x - 3;
n = 10;
d = Exponent[p0, x]; 
    a0 = Floor[Root[p0, d]];
Coeffs = 
      Table[Floor[
          Root[ k = 
              Simplify[(-x^d) (k  /.  
                      x -> (1/x + Floor[Root[k, d]]))] , 
            d]], {n}]
Frac = 
      a0 + (FromContinuedFraction[
              Coeffs])^(-1)
N[Frac]
N[Root[p0, d]]

"The test to check that adj works";

Plot[x^6 - 40  x^2 - 50, {x, -5, 5}];

Table[N[Root[x^6 - 40  x^2 - 50, i]], {i, 6}]

p0 = k = x^6 - 40  x^2 - 50;
n = 10;
d = Exponent[p0, x];
adj = 
        Catch[Do[
            If[Im[Root[p0, d - i]] == 0, Throw[i]], {i, 0, 
              d - 1}]];
a0 = Floor[Root[p0, d - adj]];
Coeffs = 
      Table[Floor[
          Root[ k = 
              Simplify[(-x^d) (k  /.  
                      x -> (1/x + Floor[Root[k, d - adj]]))] , 
            d - adj]], {n}]
Frac = 
      a0 + (FromContinuedFraction[
              Coeffs])^(-1)
N[Frac]
N[Root[p0, d - adj]]

"Notice that we don't necessarily need irreducible polys; the \
advantage may just be a shorter calculation time.";

p0 = 
        Expand[(x^6 - 40  x^2 - 50) (x + 
                2) (x + 4/5)^2];

Plot[p0, {x, -3, 5}];

k = p0;
n = 10;
d = Exponent[p0, x];
adj = 
        Catch[Do[
            If[Im[Root[p0, d - i]] == 0, Throw[i]], {i, 0, 
              d - 1}]];
a0 = Floor[Root[p0, d - adj]];
Coeffs = 
      Table[Floor[
          Root[ k = 
              Simplify[(-x^d) (k  /.  
                      x -> (1/x + Floor[Root[k, d - adj]]))] , 
            d - adj]], {n}]
Frac = 
      a0 + (FromContinuedFraction[
              Coeffs])^(-1)
N[Frac]
N[Root[p0, d - adj]]

"Nor (to find the continued fraction of the largest positive root) \
do we need just one positive root...";

p0 = 
        Expand[(x^6 - 40  x^2 - 50) (x + 2) (x + 4/5)^2  /.  
            x -> (x - 4)];

Plot[p0, {x, 0, 8}];

k = p0;
n = 10;
d = Exponent[p0, x];
adj = 
        Catch[Do[
            If[Im[Root[p0, d - i]] == 0, Throw[i]], {i, 0, 
              d - 1}]];
a0 = Floor[Root[p0, d - adj]];
Coeffs = 
      Table[Floor[
          Root[ k = 
              Simplify[(-x^d) (k  /.  
                      x -> (1/x + Floor[Root[k, d - adj]]))] , 
            d - adj]], {n}]
Frac = 
      a0 + (FromContinuedFraction[
              Coeffs])^(-1)
N[Frac]
N[Root[p0, d - adj]]

"...or even a positive leading coefficient."

p0 = 
        Expand[-(x^6 - 40  x^2 - 50) (x + 
                  2) (x + 4/5)^2  /.  
            x -> (x - 4)];

Plot[p0, {x, 0, 8}];

k = p0;
n = 10;
d = Exponent[p0, x];
adj = 
        Catch[Do[
            If[Im[Root[p0, d - i]] == 0, Throw[i]], {i, 0, 
              d - 1}]];
a0 = Floor[Root[p0, d - adj]];
Coeffs = 
      Table[Floor[
          Root[ k = 
              Simplify[(-x^d) (k  /.  
                      x -> (1/x + Floor[Root[k, d - adj]]))] , 
            d - adj]], {n}]
Frac = 
      a0 + (FromContinuedFraction[
              Coeffs])^(-1)
N[Frac]
N[Root[p0, d - adj]]

"We also don't need the positive root to be greater than one; it \
will be in all but the first poly.";

p0 = 
        Expand[x^3 - 5  x - 1  /.  x -> (x + 2)];

Plot[p0, {x, -5, 5}];

k = p0;
n = 10;
d = Exponent[p0, x];
adj = 
        Catch[Do[
            If[Im[Root[p0, d - i]] == 0, Throw[i]], {i, 0, 
              d - 1}]];
a0 = Floor[Root[p0, d - adj]];
Coeffs = 
      Table[Floor[
          Root[ k = 
              Simplify[(-x^d) (k  /.  
                      x -> (1/x + Floor[Root[k, d - adj]]))] , 
            d - adj]], {n}]
Frac = 
      a0 + (FromContinuedFraction[
              Coeffs])^(-1)
N[Frac]
N[Root[p0, d - adj]]

"Nor do we really need it to be positive, either; it will be in all \
but the first polynomial.";

p0 = 
        Expand[x^3 - 5  x - 1  /.  x -> (x + 3)];

Plot[p0, {x, -5, 5}];

k = p0;
n = 10;
d = Exponent[p0, x];
adj = 
        Catch[Do[
            If[Im[Root[p0, d - i]] == 0, Throw[i]], {i, 0, 
              d - 1}]];
a0 = Floor[Root[p0, d - adj]];
Coeffs = 
      Table[Floor[
          Root[ k = 
              Simplify[(-x^d) (k  /.  
                      x -> (1/x + Floor[Root[k, d - adj]]))] , 
            d - adj]], {n}]
Frac = 
      a0 + (FromContinuedFraction[
              Coeffs])^(-1)
N[Frac]
N[Root[p0, d - adj]]

"Finally, having the two largest roots very close to each is ok if \
the largest is separated from the other two by an integer...";

p0 = 
        Expand[x^3 - 5  x - 1   /.  
            x -> 100 (x - 1)];

Plot[p0, {x, 1, 1.2}];

k = p0;
n = 10;
d = Exponent[p0, x];
adj = 
        Catch[Do[
            If[Im[Root[p0, d - i]] == 0, Throw[i]], {i, 0, 
              d - 1}]];
a0 = Floor[Root[p0, d - adj]];
Coeffs = 
      Table[Floor[
          Root[ k = 
              Simplify[(-x^d) (k  /.  
                      x -> (1/x + Floor[Root[k, d - adj]]))] , 
            d - adj]], {n}]
Frac = 
      a0 + (FromContinuedFraction[
              Coeffs])^(-1)
N[Frac]
N[Root[p0, d - adj]]

"...but NOT ok if they have the same greatest integer";

p0 = 
        Expand[x^3 - 5  x - 1   /.  
            x -> 100 (x - 11/10)];

Plot[p0, {x, 1, 1.2}];

k = p0;
n = 10;
d = Exponent[p0, x];
adj = 
        Catch[Do[
            If[Im[Root[p0, d - i]] == 0, Throw[i]], {i, 0, 
              d - 1}]];
a0 = Floor[Root[p0, d - adj]];
Coeffs = 
      Table[Floor[
          Root[ k = 
              Simplify[(-x^d) (k  /.  
                      x -> (1/x + Floor[Root[k, d - adj]]))] , 
            d - adj]], {n}]
Frac = 
      a0 + (FromContinuedFraction[
              Coeffs])^(-1)
N[Frac]
N[Root[p0, d - adj]]

"Fortunately, by looking at the graph and adjusting the polynomial \
accordingly, such things can be easily fixed.";