As promised in a comment, here is a variant that works in fixed precision. We do get three good roots below.
First some code to do quotients of polynomials represented by their coefficient lists. This was taken from internal code for PolynomialSmithDecomposition
in some Control Theory context. (I'm allowed to do that, it was my code and written on a weekend.)
lpQuoRem[p1_, p2_] :=
Module[{p2top = p2[[-1]], top, quo, quolist, rem = p1, len, max},
top = Length[p1] - Length[p2] + 1;
If[top <= 0, Return[{{}, rem}]];
quolist = ConstantArray[0, top];
While[Length[rem] >= Length[p2], max = Max[Abs[p1]];
quo = rem[[-1]]/p2top;
quolist[[Length[rem] - Length[p2] + 1]] = quo;
rem = Most[rem] - quo*PadLeft[Most[p2], Length[rem] - 1];
len = Length[rem];
While[len > 0 && rem[[len]] == 0, len--];
rem = Take[rem, len];];
{quolist, rem}]
z[n_, c_] := If[n > 0, z[n - 1, c]^2 + c, c];
polyOrig = PolynomialQuotient[z[10, c] - z[6, c], 1 + c^2, c];
rtsPR = {};
nRoots = 3;
prec = 300;
poly = polyOrig;
Do[
Print[Timing[
Print[{Precision[poly], Exponent[poly, c]}];
aa =
FindRoot[poly, {c, I},
WorkingPrecision ->
Max[10, Min[300, Floor[Precision[poly]]]]];
fac = {c*Conjugate[c], -(c + Conjugate[c]), 1} /. aa;
Print[N@{aa, fac}];
AppendTo[rtsPR, c /. aa];
lpoly = CoefficientList[poly, c];
lpoly =
NumericalMath`FixedPrecisionEvaluate[lpQuoRem[lpoly, fac], prec];
poly = Expand[FromDigits[Reverse[lpoly[[1]]], c]];
]];
, {i, 1, nRoots}];
(* {\[Infinity],1022}
{{c->-0.000732220309309+1.00453713135 I},{1.00909538441 +0. I,0.00146444061862 +0. I,1.}}
{0.908594,Null}
{300.,1020}
{{c->-0.00956676851273+1.00673513946 I},{1.01360716408 +0. I,0.0191335370255 +0. I,1.}}
{1.250491,Null}
{300.,1018}
{{c->0.00279929953246 +1.00502058481 I},{1.01007421197 +0. I,-0.00559859906491+0. I,1.}}
{1.548567,Null} *)
Caveat: There is the possibility that eventually an accumulation of error will cause the results to be bad in the sense of not being roots to the original polynomial. A production environment would check for that and maybe raise precision when necessary. This does not require a full restart. One just polishes the "good" roots with FindRoot
to higher precision than were earlier obtained, and redoes the quotients at higher precision.
PolynomialQuotient
, which drops remainders? When I replaced it by a simple divide inside the loop, it worked fine. $\endgroup$PolynomialQuotient
is replaced by divide. Also, replacingc
byd + I
might help. Incidentally, aWorkingPrecision
of 190 gives the same results as 500, but 180 does not work. $\endgroup$$MinPrecision
in aBlock
when invokingPolynomialQuotient
. An undocumented shortcut is to useNumericalMath
FixedPrecisionEvaluate[expr,prec]`. I'll show an example where I work directly with coefficient lists in finding the quotient. $\endgroup$