Ask Your Question
1

Finding absolute values of roots of polynomials with Sage

asked 2020-12-26 20:20:00 +0200

klaaa gravatar image

updated 2020-12-27 18:46:14 +0200

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 flag offensive close merge delete

2 Answers

Sort by » oldest newest most voted
2

answered 2020-12-26 23:03:24 +0200

tmonteil gravatar image

updated 2020-12-27 00:25:42 +0200

slelievre gravatar image

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]
edit flag offensive delete link more

Comments

Why work with approximate rather than exact coefficients ?

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2020-12-27 09:46:10 +0200 )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?

klaaa gravatar imageklaaa ( 2020-12-27 18:45:01 +0200 )edit

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

Florentin Jaffredo gravatar imageFlorentin Jaffredo ( 2020-12-27 19:01:49 +0200 )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?

klaaa gravatar imageklaaa ( 2020-12-27 19:08:48 +0200 )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.

Florentin Jaffredo gravatar imageFlorentin Jaffredo ( 2020-12-27 19:22:28 +0200 )edit
1

answered 2020-12-27 20:51:03 +0200

slelievre gravatar image

updated 2020-12-27 22:57:03 +0200

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]]
edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2020-12-26 20:20:00 +0200

Seen: 76 times

Last updated: Dec 27 '20