I need to find the roots of a rational polynomial that are near i. In the following code, I try that two different ways. First, I use a constraint to only find roots in the right region. Second, I find all the roots and discard the ones not in the right region. Surprisingly, the second, wasteful method is much faster. Why?

z[0, c_] := c
z[n_, c_] := z[n - 1, c]^2 + c;
r = 1/10;
prec = 10;
expr = Expand[z[10, c] - z[6, c]];
Print[Timing[aa = c /. N[Solve[{expr == 0, Abs[c - I] < r}, c], prec];]];
Print[Timing[bb = Select[c /. N[Solve[expr == 0, c], prec], Abs[# - I] <= r &];]];

(* {4365.47,Null} {176.39,Null} *)

Clearly Solve[] is doing something terribly inefficient in its handling of constraints. I know I could use FindRoots[] to find roots one at a time and divide them out of the polynomial, then I wouldn't know when I'd found all the roots in the region and could stop looking. But when I tried that, I ran into this problem.

So how can I efficiently find all the roots within radius r of i?

And what the heck is going with Solve[]?

  • $\begingroup$ Check out CountRoots. It can count roots in a square region of the complex plane. $\endgroup$
    – John McGee
    Commented Jan 25, 2016 at 11:25
  • $\begingroup$ Why not directly use SolveValues[…, c] (instead of c /. Solve[…])? $\endgroup$
    – user688486
    Commented Dec 25, 2023 at 9:54
  • $\begingroup$ @user688486 Because I've never heard of that function because Mathematica documentation is terrible. $\endgroup$ Commented Dec 26, 2023 at 13:59

2 Answers 2


Here is an entirely different numeric method, from this Usenet post. The idea is to use the Cauchy Integral formula for root counting, isolation, and refinement. I think the original code had some problems so below is modified (and thus probably has different problems..)

SetOptions[NIntegrate, WorkingPrecision -> 150, AccuracyGoal -> 8, 
  PrecisionGoal -> 8];

eps = 1.3*10^(-6);
isolateRootsInRectangle[func_, a_, b_, c_, d_, z_, eps_] := 
 Module[{numroots, avgx = (a + b)/2, avgy = (c + d)/2, r1, r2, r3, r4,
    len = 0, allrts}, 
  numroots = 
   Round[rectangleIntegrate[func, a - eps, b + eps, c - eps, d + eps, 
  If[numroots <= 0, Return[{}]];
  If[numroots == 1, Return[{{a, b, c, d}}]];
  r1 = isolateRootsInRectangle[func, a, avgx, c, avgy, z, eps];
  r2 = isolateRootsInRectangle[func, a, avgx, avgy, d, z, eps];
  r3 = isolateRootsInRectangle[func, avgx, b, c, avgy, z, eps];
  r4 = isolateRootsInRectangle[func, avgx, b, avgy, d, z, eps];
  Join[r1, r2, r3, r4]]

I now find intervals for roots that lie in a box of half-length 1/10 centered at I.

r = 1/10;
Timing[intervals = 
   isolateRootsInRectangle[expr, -r, r, 1 - r, 1 + r, c, eps];]

(* Out[567]= {27.433912, Null} *)

This may have duplicates due to use of slight overlaps on region, and I did not attempt to remove them.

Now we refine them. Could use FindRoot with start points at interval centers, but this is more reliable (and a FindRoot based on these values might be better still).

Timing[rts = refineRoots[expr, intervals, c, eps];]

(* Out[573]= {6.487449, Null} *)

As for working with a circular rather than square region, there is more code at the link above. But that involves extra fuss and it's probably simpler just to remove those that are inside the box but outside the circle in question.


Solve also finds all complex roots and selects those that satisfy the inequality condition. Most of the timing difference is due to Solve evaluating the condition with exact values of roots and your code evaluating the condition with 10-digit approximations of roots. Exact evaluation of inequalities computes the expressions involved to at least 54 bits of precision. In this example computing the extra digits of precison is costly.

In[3]:= (rts=Solve[expr == 0, c]);//Timing
Out[3]= {126.159, Null}

In[4]:= (sol=Select[c /. rts, Abs[# - I] <= r &]);//Timing
Out[4]= {1229.03, Null}
  • $\begingroup$ Aha, that is a helpful insight on the why part. Thank you. Any thoughts on the "So how can I efficiently find all the roots within radius r of i?" part? $\endgroup$ Commented Jan 26, 2016 at 2:20
  • $\begingroup$ I think your method might be optimal for this polynomial. Most of the Solve time is spent factoring the polynomial. In[12]:= (fl=FactorList[expr]);//Timing Out[12]= {112.266, Null} $\endgroup$ Commented Jan 26, 2016 at 17:32
  • 1
    $\begingroup$ The factors have degrees <= 384. I don't think numerically finding roots of the original polynomial of degree 1024 would be faster. $\endgroup$ Commented Jan 26, 2016 at 17:41
  • $\begingroup$ I don't think this could be optimal, because it's wasting find find ALL the roots instead of stopping when the roots of interest have all be found $\endgroup$ Commented Jan 26, 2016 at 19:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.