Finding absolute values of roots of polynomials with Sage

polynomials

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

klaaa
671 ●3 ●14 ●38

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

add a comment

1

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

slelievre
17654 ●22 ●160 ●348 http://carva.org/samue...

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

add a comment

2

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

tmonteil
27323 ●31 ●202 ●514 http://wiki.sagemath.o...

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

slelievre
17654 ●22 ●160 ●348 http://carva.org/samue...

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

Comments

Why work with approximate rather than exact coefficients ?

Emmanuel 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 ( 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 ( 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 ( 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 ( 2020-12-27 19:22:28 +0200 )edit

add a comment

Finding absolute values of roots of polynomials with Sage

2 Answers

Roots of monic integer polynomials

Polynomials with complex coefficients

Monic integer polynomials

Comments

Your Answer

Question Tools

Stats

Related questions

Finding absolute values of roots of polynomials with Sage edit

2 Answers

Roots of monic integer polynomials

Polynomials with complex coefficients

Monic integer polynomials

Comments

Your Answer

Question Tools

Stats

Related questions

Finding absolute values of roots of polynomials with Sage