I need to find all the rational roots of a bunch of high-order polynomials with rational complex coefficients. Is there an efficient way to do this?

My best effort on an example polynomial:

expr = (12 x - 29)*(31 x + 113)*(501 x - Prime[128])*(x^4 - 7 x^2 + 1001 x - 20)*
       (1001*x^5 + 47*x^4 + x^3 + 91*x^2 + 144 x - 1001) // Expand;

If[IrreduciblePolynomialQ[expr], r = {},
  a = Factor[expr];
  n = Count[Table[Exponent[a[[i]], x], {i, 1, Length[a]}], 1];
  r = Table[Solve[a[[i]] == 0, x][[1]], {i, 1, n}]];

(*  {{x -> 29/12}, {x -> -(113/31)}, {x -> 719/501}}  *)

The problem with this method are: 1) It does a full Factor[] operation on the polynomial when there's only a specific type of factor I care about 2) In light of this discussion I am not certain that Factor[] will successfully find all 1st-order polynomial factors.

  • 2
    $\begingroup$ Your example might not include "bad cases", perhaps because it doesn't include either rational or complex coefficients. BUT if ignore that and apply the rational root theorem, by constructing all integer factors of 47170383260 and of 186558372 and forming all rationals from all subsets of those two lists then I almost instantly find exactly and only the rational roots of your expr. Does this then apply in general to your problem?!? $\endgroup$
    – Bill
    Commented May 3, 2016 at 4:27
  • $\begingroup$ @Bill My example does include rational coefficients, so I'm not sure what you meant by that. And I thought of the rational roots approach you suggest, but it gets so complicated for complex coefficients that I started again to hope there was a more efficient way. And even if that is the most efficient way, I would wonder about how to code it efficiently. $\endgroup$ Commented May 3, 2016 at 5:13
  • $\begingroup$ How exactly do you suppose Factor will miss a linear factor? I might give the general question some thought, except it seems to be based on an utterly bizarre premise. $\endgroup$ Commented May 3, 2016 at 15:22
  • $\begingroup$ @DanielLichtblau Well first, I'm 99.999% sure it won't miss a linear factor, but I don't like to assume. Second, even if my code above works 100% of the time, my question is still: can this calculation be done more efficiently? I've learned a lot about MMa and good coding practice by asking questions like that. $\endgroup$ Commented May 3, 2016 at 16:54
  • $\begingroup$ Factor works correctly to identify all factors over the rationals. This includes of course those that happen to be linear. As for possibly more efficient methods, yes, they exist. I might post one later if time allows. $\endgroup$ Commented May 3, 2016 at 16:57

1 Answer 1


I'm going to restrict to the case of rational coefficients. There are ways to extend to complex rational coefficients but that's more than I have time or desire to consider right now.

I'll illustrate an efficient methodology with this example. Along the way I will say a bit about modest improvements that can be made. We'll start with the polynomial and the desired result.

expr = 
 Expand[(12 x - 29)*(31 x + 113)*(501 x - Prime[128])*(x^4 - 7 x^2 + 
     1001 x - 20)*(1001*x^5 + 47*x^4 + x^3 + 91*x^2 + 144 x - 1001)]
ratrts = Cases[x /. Solve[expr == 0, x], _Rational]

(* Out[836]= 47170383260 - 2407110038803 x + 2331054513893 x^2 - 
 39602581260 x^3 - 374204769604 x^4 + 91472770391 x^5 + 
 2323640009612 x^6 - 1991681044106 x^7 - 19601204064 x^8 + 
 189236294145 x^9 - 3279862609 x^10 - 29790027 x^11 + 186558372 x^12

Out[837]= {-(113/31), 719/501, 29/12} *)

Step 1: Choose a prime pr such that the leading coefficient does not vanish modulo pr and moreover there are no multiple roots. I omit the check for those features.

Step 2: Find a bound on size of numerator and denominator of any rational root. By the factor theorem of basic algebra this is just the first and last coefficients of the polynomial (assuming last is not zero, since any integer divides zero).

Step 3: Determine the first integral power of pr to exceed twice that bound squared (I think it would suffice here just to use twice the product of the coefficients).

pr = Prime[123457];
trcoeff = Coefficient[expr, x, 0];
ldcoeff = Coefficient[expr, x, Exponent[expr, x]];
liftheight = 2*Max[ldcoeff, trcoeff]^2;
liftpower = 1;
ppow = pr^liftpower;
While[ppow < liftheight,
  liftpower *= 2; ppow = ppow^2];

Step 4: Take GCD of the polynomial with x^pr-x mod pr. This removes any nonlinear factors modulo pr.

Step 5: Extract the factors of what is left. This could be done with Berlekamp's algorithm or otherwise (I think the main alternative is due to Cantor). But could simply use Roots or even Factor since they will end up using one of those under the hood with relatively little overhead.

p2 = PolynomialGCD[PolynomialMod[x^pr, {expr, pr}] - x, expr, 
   Modulus -> pr];
prroots = x /. {ToRules[Roots[p2 == 0, x, Modulus -> pr]]}

(* Out[853]= {400212, 952535, 998130, 1000814, 1116416, 1133236, 1352614} *)

Step 6: Lift each root using Hensel's lemma. That is to say, for each root correct modulo pr, we find a correction term so that it becomes valid mod pr^2. Iterate that squaring/correcting up to the bound found in step 3 since that's what is needed to do rational recovery.

Step 7: Recover rational values from power-of-prime lifted factors.

Step 8; Remove all results that exceed the size of valid rational roots.

We combine these steps below. The code for lift I'm being lazy and using Solve but it's easy to use e.g. CoefficientList and related in order to get the corrected roots at each lifting step.

lift[ee_, x_, x0_, p_] := Module[
  {deriv = D[ee, x], dx, corr},
  deriv = D[ee, x] /. x -> x0;
  corr = Expand[PolynomialMod[(ee + deriv*p*dx) /. x -> x0, p^2]/p];
  corr = dx /. First[Solve[corr == 0, dx, Modulus -> p]];
  x0 + corr*p

rationalRecover[x_, pk_] :=
 ((#[[2, 2]]/#[[1, 2, 2]]) &)[Internal`HGCD[pk, x]]

iftedrts = Table[
   liftrt = rt;
   ppow = pr;
   Do[liftrt = lift[expr, x, liftrt, ppow]; 
    ppow = ppow^2, {j, liftpower - 1}];
   rationalRecover[liftrt, ppow]
   , {rt, prroots}];
 Abs[Numerator[#]] <= trcoeff && Abs[Denominator[#]] <= ldcoeff &]

(* Out[834]= {719/501, -(113/31), 29/12} *)

For small examples like this I don't see any advantage over just using Factor. At large degree, with large coefficients, and perhaps but few rational roots, it would be a different matter entirely.


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.