How to efficiently find only the rational roots of a rational complex polynomial?

Question

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}]];
r

(*  {{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.

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}];
Select[liftedrts, 
 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.