Finding absolute values of roots of polynomials with Sage

Given a monic integer polynomial $p(x)$ in one variable $x$. Is it possible using Sage to obtain all absolute values of the roots of $p(x)$ in a list?

(Or at least all absolute values of roots that have absolute value (say) at most 2, if there general case is not possible and one needs a (small) bound.)

In my work there appear polynomials with large degree up to 40-50 so Im not sure whether such a thing is quickly possible for nowadays computers. But a solution would also be interesting even when it works only for degrees up to 20 or so.

The solution should give the correct answer for the absolute value of the roots for 5 or 6 decimal places.

edit retag close merge delete

Sort by ยป oldest newest most voted

Polynomials with complex coefficients

You can define the ring of (numerical) polynomials over C as follows:

sage: R.<x> = CDF[]
sage: R
Univariate Polynomial Ring in x over Complex Double Field


You can construct a random polynomial of degree 30 as follows:

sage: P = R.random_element(degree=30)
sage: P
(0.7368648392813568 + 0.7644659236040987*I)*x^30 + (-0.9521231496916671 - 0.5668264771579159*I)*x^29 + (-0.5602325090820397 + 0.08373185660033844*I)*x^28 + (0.11778992480095218 - 0.2031034081348333*I)*x^27 + (0.3516830229409873 + 0.2968075058912192*I)*x^26 + (-0.7794485386929075 + 0.35282724454074765*I)*x^25 + (0.9529334155548153 - 0.10858391867189021*I)*x^24 + (0.4346153404358597 - 0.7919849906353298*I)*x^23 + (-0.4122547991033889 - 0.09369196690243964*I)*x^22 + (0.9072745543379173 - 0.93084223388219*I)*x^21 + (0.0917339388366103 + 0.8400827723548798*I)*x^20 + (0.620441590530908 + 0.5497959418695815*I)*x^19 + (0.08961053619498971 - 0.3245694231244378*I)*x^18 + (-0.1573830771827882 - 0.0999821789852311*I)*x^17 + (-0.42868409209771885 + 0.17826524252278442*I)*x^16 + (0.08423322767931607 - 0.7883353014539329*I)*x^15 + (-0.4854333236724999 + 0.14977131931558807*I)*x^14 + (-0.44684046421460066 - 0.46310410665013624*I)*x^13 + (0.6526594274639133 - 0.4196336531832292*I)*x^12 + (-0.27787355156700233 + 0.10848474232286232*I)*x^11 + (-0.11101383721798008 - 0.5037475397998767*I)*x^10 + (-0.13186950326025015 - 0.4838574269429168*I)*x^9 + (0.8187926299228285 + 0.8711795580308974*I)*x^8 + (-0.6860149090655201 - 0.9849078620505949*I)*x^7 + (-0.6918369434398146 - 0.9814493553853534*I)*x^6 + (0.3499220158155265 + 0.5747414242705378*I)*x^5 + (-0.893601061756806 + 0.5607846129706036*I)*x^4 + (-0.12889260181776785 + 0.7202322287337615*I)*x^3 + (-0.6940468382565192 - 0.5557200414145016*I)*x^2 + (-0.08633744827278456 + 0.7424063585043299*I)*x + 0.6230123145504038 + 0.6286232382291024*I


You can find its roots as follows:

sage: P.roots(multiplicities=False)
[-1.0120138172468134 + 0.20426669296859173*I,
-0.962581077360827 - 0.5296864154447034*I,
-0.9364728333831983 - 0.24073599552387334*I,
-0.9201821425300504 + 0.002487918823541706*I,
-0.8599059855418224 + 0.492490723674138*I,
-0.8215080191130956 - 0.6427711927416188*I,
-0.7327655570630038 + 0.6808285540382258*I,
-0.7017713764174295 + 0.9067735719613969*I,
-0.6414579210848824 + 0.1305690899208492*I,
-0.5346277142391518 - 0.7813046856018692*I,
-0.36837899733455526 - 0.896830748646684*I,
-0.3534478170558439 + 0.7077101949275129*I,
-0.28232967917268137 + 0.9677582197606973*I,
-0.22178572802727103 - 0.9734268214408631*I,
0.015053958533473745 + 1.035138545876896*I,
0.09170010634309536 - 1.0719575404124397*I,
0.23422262713071448 + 0.9812411153412857*I,
0.243546812209695 - 1.029432097206379*I,
0.39934906552852045 - 0.8111814833649242*I,
0.4716233302309813 + 0.8546354765171169*I,
0.582436456639921 - 0.7624604937976489*I,
0.6491856986029882 + 0.6138112252444622*I,
0.7124110338205422 + 0.4351851664034465*I,
0.7570474447707015 + 0.7787758630291586*I,
0.831513193382847 - 0.6444495428445836*I,
0.9648971684027969 - 0.38704090145535897*I,
0.9660222981279798 + 0.24237600403657947*I,
0.9892771693661305 + 0.06311132544291641*I,
1.1525531039685055 + 0.017246602835138347*I,
1.2950663367376647 - 0.618272306837997*I]


You can construct the list of its absolute values as follows:

sage: [abs(r) for r in P.roots(multiplicities=False)]
[1.0324228049373916,
1.0986946933519746,
0.9669204658117574,
0.9201855058471996,
0.9909517732340606,
1.0430868763845154,
1.0002364128602976,
1.1466130016559397,
0.6546117565132008,
0.946712102265528,
0.9695403433562648,
0.7910619946539459,
1.0081002031822746,
0.9983730193955808,
1.0352480045031056,
1.0758726123247224,
1.0088083888906871,
1.0578494658958941,
0.9041543425171509,
0.9761303003901686,
0.9594676808655002,
0.8934239147833227,
0.8348147160698889,
1.0860997553036549,
1.052012074095318,
1.0396283975493898,
0.995964461119574,
0.9912882311559164,
1.152682134318386,
1.4350809949107028]


Monic integer polynomials

Likewise define the ring of polynomials in x with integer coefficients:

sage: Zx = ZZ['x']
sage: x = Zx.gen()


Pick a random monic polynomial of degree 56:

sage: q = Zx.random_element(degree=56)  # not necessarily monic
sage: p = (1 - q[56])*x^56 + q  # now monic by making leading coefficient 1
x^56 - x^55 + x^54 - 7*x^53 - 2*x^52 - x^51 + 5*x^50 - x^49 - x^48 + x^47
- 3*x^44 + 3*x^43 + 3*x^42 + 20*x^41 - x^40 + x^39 + 125*x^38 + 9*x^37
- x^36 - 2*x^35 + 3*x^34 + 10*x^32 + x^30 + 99*x^29 + 7*x^27 + 13*x^25
+ x^23 + x^22 + 2*x^21 + 3*x^20 + 2*x^19 + x^17 + 2*x^14 - x^12 + 5*x^11
+ x^10 - x^8 + 7*x^7 - x^6 + x^5 - 5*x^4 - x^3 + 5*x^2 + x - 1


Get the sorted list of absolute values of its complex roots:

sage: sorted(set(abs(z) for z in p.complex_roots()))
[0.392775879097575, 0.618313001041181, 0.829872478969433,
0.876156325725046, 0.877139969220623, 0.896140603140608,
0.897355058131500, 0.898822808765931, 0.902218110575174,
0.908268048474002, 0.909765827992206, 0.914374431830909,
0.923321003752644, 0.935403911031256, 0.938567671569716,
0.970873370350356, 0.972725264654593, 0.975123271228316,
0.977132068373407, 1.01049931640877, 1.15177593823796,
1.16627154233445, 1.18058797182489, 1.21320658003016,
1.21544302417961, 1.23985793698233, 1.26758074663977,
1.28869109993426, 1.80658377491657, 2.17180464952677]

more

Why work with approximate rather than exact coefficients ?

( 2020-12-27 02:46:10 -0600 )edit

Thank you but this is not very precise. For example:

R.<x> = CDF[]
P=x^7 + 7*x^6 + 21*x^5 + 35*x^4 + 35*x^3 + 21*x^2 + 7*x + 1
UU=P.roots()
display(UU)


this polynomial has only the zero -1 (it is equal to (x+1)^7). But the approximated zeros are too imprecise. Is there a better way to get zeros that are correct at least until the 5. or 6. decimal place?

( 2020-12-27 11:45:01 -0600 )edit

You can always replace CDF[] by PolynomialRing(ComplexField(100)) to work with 100 bits of precision. In my case anything above 6 works.

( 2020-12-27 12:01:49 -0600 )edit

Thank you very much! Do you have a realistic expectation what to use to get a fast and correct result that is correct up to 7 decimals? For example is PolynomialRing(ComplexField(100)) good enough for that?

( 2020-12-27 12:08:48 -0600 )edit

It probably depends on the polynomial. My guess is that it will become worse with increasing degree. One way to find the appropriate precision is to compare the result with a known root. Here for example I increase the precision in order to compute sqrt(2) with 7 digits of precision.

for i in range(1, 1000):
R.<x> = PolynomialRing(ComplexField(i))
P = R((x-sqrt(2))^7)
UU = P.roots()
if all([abs(x[0]-sqrt(2)).n() < 1e-7 for x in UU]):
print(i)
break


172 in this case.

( 2020-12-27 12:22:28 -0600 )edit

Roots of monic integer polynomials

Define the ring of polynomials in x with integer coefficients:

sage: Zx = ZZ['x']
sage: x = Zx.gen()


Pick a random monic polynomial of degree 56:

sage: q = Zx.random_element(degree=56)  # not necessarily monic
sage: p = (1 - q[56])*x^56 + q  # now monic by making leading coefficient 1
sage: p
x^56 + 27*x^54 + x^53 + 2*x^52 - 4*x^51 + 6*x^49 - 3*x^48 - 2*x^47
+ 3*x^46 - 3*x^45 + x^43 + x^41 - x^38 + x^36 - 3*x^34 - 2*x^33
+ 14*x^32 - 2*x^31 - x^30 + x^29 + x^28 + x^27 - x^26 - 3*x^25
- 4*x^24 + x^23 - 13*x^20 + x^19 - 3*x^18 - 3*x^16 + x^15 - 21*x^13
+ x^11 + x^10 - 2*x^8 - 2*x^6 + 4*x^5 + 6*x^4 - x^3 + 3*x^2 + x - 2


Compute the roots as complex balls:

sage: rr = p.roots(CBF, multiplicities=False)
sage: len(rr)
56


Take absolute values, excluding those with negative imaginary part so as to take only one per pair of conjugate non-real roots, and sort:

sage: arr = sorted(abs(r) for r in rr if r.imag() >= 0)
sage: len(arr)
26
sage: arr
[[0.790151256486441 +/- 2.55e-16], [0.834175531654164 +/- 4.27e-16],
[0.867259115927234 +/- 5.27e-16], [0.938405173485585 +/- 2.82e-16],
[0.951240193645357 +/- 5.04e-16], [0.959393543021787 +/- 4.84e-16],
[0.968049174961464 +/- 3.74e-16], [0.970884838146261 +/- 3.72e-16],
[0.971668448923249 +/- 5.60e-16], [0.981855232640659 +/- 3.31e-16],
[0.988836704194625 +/- 3.35e-16], [0.989893894583637 +/- 3.10e-16],
[0.992900145032832 +/- 2.58e-16], [0.993485276689509 +/- 4.86e-16],
[0.995962550963240 +/- 5.19e-16], [0.998272030405643 +/- 5.03e-16],
[0.998888797451878 +/- 3.15e-16], [1.001077438291740 +/- 2.41e-16],
[1.003483438283551 +/- 3.24e-16], [1.003920944016440 +/- 5.09e-16],
[1.005628887116308 +/- 5.22e-16], [1.011338025074980 +/- 5.16e-16],
[1.011575459414411 +/- 3.84e-16], [1.012926800214650 +/- 5.86e-16],
[1.035073067121418 +/- 3.44e-16], [5.18922034049960 +/- 2.92e-15]]

more