I think you can avoid the Newton-Raphson altogether, since cubic is solvable in radicals.
Here is the complete algorithm (working under constraints outlined in the problem), in Mathematica:
Clear[PositiveCubicRoot];
PositiveCubicRoot[p_, q_, r_] :=
Module[{po3 = p/3, a, b, det, abs, arg},
b = ( po3^3 - po3 q/2 + r/2);
a = (-po3^2 + q/3);
det = a^3 + b^2;
If[det >= 0,
det = Power[Sqrt[det] - b, 1/3];
-po3 - a/det + det
,
(* evaluate real part, imaginary parts cancel anyway *)
abs = Sqrt[-a^3];
arg = ArcCos[-b/abs];
abs = Power[abs, 1/3];
abs = (abs - a/abs);
arg = -po3 + abs*Cos[arg/3]
]
]
Then
In[222]:= PositiveCubicRoot[-2.52111798, -71.424692, -129.51520]
Out[222]= 10.499
However, if the Newton-Raphson method must be used, then a good initial guess is imperative. Binary division is a good method to isolate the root in the case at hand.
We start with an arbitrary point $x$, I chose $x=1$, and double it while the polynomial
at that point is negative. Then do binary division a certain number of times (the code below does it twice). Ultimately, polish it off with Newton-Raphson iterations:
In[283]:=
NewtonRaphsonStartingPoint[{p_, q_, r_}] := Module[{x1=0, x2=1,f1,f2,xm,fm},
f1 = r + x1 (q + (p + x1) x1);
While[(f2 = r + x2 (q + (p + x2) x2)) <= 0,
x1 = x2; f1 = f2; x2 = 2 x2];
Do[xm = (x1 + x2)/2; fm = r + xm (q + (p + xm) xm);
If[fm <= 0, f1 = fm; x1 = xm, f2 = fm; x2 = xm], {i, 2}];
(f2 x2 - f1 x1)/(f2 - f1)
];
NewtonRaphsonIterate[{p_, q_, r_}, x0_Real] :=
FixedPoint[
Function[x, x - (r + x (q + (p + x) x))/(q + (2 p + 3 x) x)], x0]
In[285]:=
NewtonRaphson[p_, q_, r_] :=
NewtonRaphsonIterate[{p, q, r}, NewtonRaphsonStartingPoint[{p, q, r}]]
In[286]:= NewtonRaphson[-2.52111798, -71.424692, -129.51520]
Out[286]= 10.499
There is such a formula: consider
$$\frac{x_n+\sqrt y}{x_n-\sqrt y}=\frac{\frac{x_{n-1}^2+y}{2x_{n-1}}+\sqrt y}{\frac{x_{n-1}^2+y}{2x_{n-1}}-\sqrt y}=\frac{(x_{n-1}+\sqrt y)^2}{(x_{n-1}-\sqrt y)^2}=\left(\frac{x_{n-1}+\sqrt y}{x_{n-1}-\sqrt y}\right)^2.$$
By recurrence,
$$\frac{x_n+\sqrt y}{x_n-\sqrt y}=\left(\frac{x_{0}+\sqrt y}{x_{0}-\sqrt y}\right)^{2^n}.$$
If you want to achieve $2^{-b}$ relative accuracy, $x_n=(1+2^{-b})\sqrt y$,
$$2^n=\frac{\log_2\frac{(1+2^{-b})\sqrt y+\sqrt y}{(1+2^{-b})\sqrt y-\sqrt y}}{\log_2\left|\frac{x_{0}+\sqrt y}{x_{0}-\sqrt y}\right|},$$
$$n=\log_2\left(\log_2\frac{2+2^{-b}}{2^{-b}}\right)-\log_2\left(\log_2\left|\frac{x_{0}+\sqrt y}{x_{0}-\sqrt y}\right|\right).$$
The first term relates to the desired accuracy. The second is a penalty you pay for providing an inaccurate initial estimate.
If the floating-point representation of $y$ is available, a very good starting approximation is obtained by setting the mantissa to $1$ and halving the exponent (with rounding). This results in an estimate which is at worse a factor $\sqrt 2$ away from the true square root.
$$n=\log_2\left(\log_2\left(2^{b+1}+1\right)-\log_2\left(\log_2\frac{\sqrt 2+1}{\sqrt 2-1}\right)\right)
\approx\log_2(b+1)-1.35.$$
In the case of single precision (23 bits mantissa), 4 iterations are always enough. For double precision (52 bits), 5 iterations.
On the opposite, if $1$ is used as a start and $y$ is much larger, $\log_2\left|\frac{1+\sqrt y}{1-\sqrt y}\right|$ is close to $\frac{2}{\ln(2)\sqrt y}$ and the formula degenerates to
$$n\approx\log_2(b+1)+\log_2(\sqrt y)-1.53.$$
Quadratic convergence is lost as the second term is linear in the exponent of $y$.
Best Answer
If you use Halley's method, you exhibit cubic convergence! This method is second in the class of Householder's methods.
Halley's method is: $$ x_{n+1} = x_n - \frac{2f(x_n)f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)} $$ If we let $$f(x) = x^2 - a$$ which meets the criteria, (continuous second derivative)
Then Halley's method is:
$$ x_{n+1} = x_n - \frac{\left(2x_n^3 - 2ax_n\right)}{3x_n^2 + a} $$ Which has the simplification: $$ x_{n+1} = \frac{x_n^3 + 3ax_n}{3x_n^2 + a} $$ I also will add this document which discusses extensions of the newtonian method.
There exists an extension due to Potra and Pták called the “two-step method” that may be re-written as the iterative scheme $$x_{n+1} =x_n − \frac{f(x_n)+f\left(x_n − \frac{f(x_n)}{f′(x_n)}\right)}{f'(x_n)}$$ that converges cubically in some neighborhood of of the root $x^{*}$ which does not require the computation of the second derivative.
See: On Newton-type methods with cubic convergence for more information on this topic.
As Hurkyl and others have noted, your best bet is to just use Newton's Method. These alternative methods generally come with more operations per iteration. They aren't really worth the computational cost, but they are a good comparison.