The following equation solving too slow (compared with Wolfram Mathematica)

[stevenleeS0ht]

The following expression:

import sympy

sympy.solve(sympy.Eq(sympy.S('59.0*(1 - 1/(ytm/2 + 1)**16)/(sqrt((ytm/2 + 1)**2) - 1) + 1000/(ytm/2 + 1)**16'), 779.92))

Take so long time, I cannot wait for it to return result (after several hours, it still didn't return result).

If I try to solve it in Mathematica Solve[1000/(1+ytm/2)^16+(118 (1-1/(1+ytm/2)^16))/ytm == 779.92], the result will be returned before I rise my Return key.

[oscarbenjamin]

It's better to use nsolve for something like this:

In [108]: ytm = Symbol('ytm')

In [109]: eq = Eq(59.0*(1 - 1/(ytm/2 + 1)**16)/(sqrt((ytm/2 + 1)**2) - 1) + 1000/(ytm/2 + 1)**16, 779.92)

In [110]: nsolve(eq, ytm, 0.1)
Out[110]: 0.169199699318998

In [112]: %time nsolve(eq, ytm, 0.1)
CPU times: user 12.5 ms, sys: 89 µs, total: 12.6 ms
Wall time: 12.8 ms
Out[112]: 0.169199699318998

What kind of output does Mathematica return?

The reason this is slow is because after unrad the solutions are found as the roots of this poly:

Poly(95043001*ytm**65 + 12545676132*ytm**64 + 815029081851*ytm**63 + 34737038864768*ytm**62 + 1092443181588864*ytm**61 + 27033945759733248*ytm**60 + 548199666379615488*ytm**59 + 9367045534245648384*ytm**58 + 137636261949947277312*ytm**57 + 1766205976510378213376*ytm**56 + 20035065680637627279360*ytm**55 + 202865675356690523455488*ytm**54 + 1848230840737496175837184*ytm**53 + 15251476819261937637261312*ytm**52 + 114630160156436181653323776*ytm**51 + 788457483403739474714689536*ytm**50 + 4983305790039208809504440320*ytm**49 + 29043065170208026970643496960*ytm**48 + 156558597284899313010955386880*ytm**47 + 782661702029057497437348823040*ytm**46 + 3636979493736225855303330037760*ytm**45 + 15741914847465179713116978544640*ytm**44 + 63576080853266087524914293637120*ytm**43 + 239950153855567263056140393512960*ytm**42 + 847468399712575035118008643092480*ytm**41 + 2804183176035662026628909366247424*ytm**40 + 8701675686667671954144473443729408*ytm**39 + 25344288757181483747782371829612544*ytm**38 + 69333560162977193120983178377428992*ytm**37 + 178254569984163413657835278031650816*ytm**36 + 430887839096586051718601028866670592*ytm**35 + 979611747325145040689863847451820032*ytm**34 + 2095069232453211474522379340927205376*ytm**33 + 4215367358157470666141580126026989568*ytm**32 + 7979082183629421417575423514681278464*ytm**31 + 14206534464374020218282232758633758720*ytm**30 + 23786290591698389089330153768607547392*ytm**29 + 37436997193116426159405230419771129856*ytm**28 + 55358911427031235467068320075122475008*ytm**27 + 76860567715705247183097868601350160384*ytm**26 + 100115791956243635409507646298394722304*ytm**25 + 122227308519736044769153215122940887040*ytm**24 + 139703390290878479441799343066180485120*ytm**23 + 149292487527700125374497542678889103360*ytm**22 + 148929558810059501001683580071318650880*ytm**21 + 138432976101339414393432658848902021120*ytm**20 + 119641184242627099279216777800586362880*ytm**19 + 95897303922420257858095850789577687040*ytm**18 + 71075012380718275930817278610054840320*ytm**17 + 48535635480348516776384968068861788160*ytm**16 + 30406373367009521823221596708556242944*ytm**15 + 17383407659186303251849539610988249088*ytm**14 + 9009705731771310824491496571177795584*ytm**13 + 4197939935555399399050819212633178112*ytm**12 + 1738949892589644498724259893950283776*ytm**11 + 630664358453423139656279407766536192*ytm**10 + 195758173766356235005063816818982912*ytm**9 + 50101117415821219039109970307055616*ytm**8 + 9817722026961490489485534154457088*ytm**7 + 1183752679154209281091930658177024*ytm**6 - 27513872814296384523330387968000*ytm**5 - 50422798089985294691191245242368*ytm**4 - 11688205059328477916116627226624*ytm**3 - 1307430557707125370964810924032*ytm**2 - 45757297558507961693279092736*ytm + 2953677977474090731950833664, ytm, domain='ZZ')

The poly.all_roots method can easily factor this but it then tries to create CRootOf for this factor

PurePoly(9749*x**16 + 310493*x**15 + 4632320*x**14 + 42967520*x**13 + 277282880*x**12 + 1319724224*x**11 + 4790309888*x**10 + 13519700480*x**9 + 29960353280*x**8 + 52242910720*x**7 + 71304183808*x**6 + 75115995136*x**5 + 59481210880*x**4 + 33728020480*x**3 + 12400721920*x**2 + 2211315712*x - 953614336, x, domain='ZZ')

which is then slow because the root isolation code is slow for complex roots. I don't understand why the root isolation is so slow because I know that it is possible to make something a lot faster.

It would be better to use Poly.nroots but internally solve with rational=True (the default) does not know that the user originally provided floats:

In [8]: %time PurePoly(9749*x**16 + 310493*x**15 + 4632320*x**14 + 42967520*x**13 + 277282880*x**12 + 1319724224*x**11 + 4790309888*x**10 + 13519700480*x**9 +
   ...:  29960353280*x**8 + 52242910720*x**7 + 71304183808*x**6 + 75115995136*x**5 + 59481210880*x**4 + 33728020480*x**3 + 12400721920*x**2 + 2211315712*x - 9
   ...: 53614336, x, domain='ZZ').nroots()
CPU times: user 225 ms, sys: 1.33 ms, total: 226 ms
Wall time: 228 ms
Out[8]: 
[-4.03029835203001, 0.169199699318998, -3.87567487140385 - 0.777142968354318⋅ⅈ, -3.87567487140385 + 0.777142968354318⋅ⅈ, -3.43534495188389 - 1.4359244411791
2⋅ⅈ, -3.43534495188389 + 1.43592444117912⋅ⅈ, -2.77634665334044 - 1.8759936941525⋅ⅈ, -2.77634665334044 + 1.8759936941525⋅ⅈ, -1.99901107918184 - 2.03027332569
222⋅ⅈ, -1.99901107918184 + 2.03027332569222⋅ⅈ, -1.22169322936442 - 1.87513320184366⋅ⅈ, -1.22169322936442 + 1.87513320184366⋅ⅈ, -0.562780544415769 - 1.433849
65693368⋅ⅈ, -0.562780544415769 + 1.43384965693368⋅ⅈ, -0.122950559563576 - 0.772171826498919⋅ⅈ, -0.122950559563576 + 0.772171826498919⋅ⅈ]

[oscarbenjamin]

This kind of equation could be solved much faster if solve only tried to find real solutions. Many equations like this reduce to large polynomials and then finding all real roots of large polynomials (with RootOf) is much faster than trying to isolate all of the complex roots.