3
$\begingroup$

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.

$\endgroup$
  • 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 May 3 '16 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$ – Jerry Guern May 3 '16 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$ – Daniel Lichtblau May 3 '16 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$ – Jerry Guern May 3 '16 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$ – Daniel Lichtblau May 3 '16 at 16:57
3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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