HP Forum Archive 20

[ Return to Index | Top of Index ]

Any expert for polynomial roots here? Message #1 Posted by fhub on 8 Nov 2011, 7:40 a.m. For my WP34s Polynomial RootSolver I've used the Bairstow method, but I also checked many other algorithms that I've found. One of these is the Laguerre method - here's a short WikiPedia article about it: http://en.wikipedia.org/wiki/Laguerre%27s_method It's stated there that this Laguerre method should be VERY powerful in 2 ways: it converges almost always (independent of the initial guesses), and it converges very quickly (usually 3rd order convergence). Now I made a 2nd version of my PolyRoots where I replaced the Bairstow method by this Laguerre method, and I'm absolutely disappointed: although it works correctly (i.e. finds the correct roots), it is damned slow - in average it needs about 3-5 times as many iterations as the Bairstow method (although it SHOULD be much better)! :-( Has anyone here experiences with polynomial rootsolving methods? Maybe someone could tell me what might be wrong with this Laguerre method? (I'm talking about polynomials with real coefficients - maybe this is the reason why Laguerre works so badly?). Franz Edited: 8 Nov 2011, 7:58 a.m. after one or more responses were posted

Re: Any expert for polynomial roots here? Message #2 Posted by Valentin Albillo on 8 Nov 2011, 7:52 a.m., in response to message #1 by fhub In case you haven't already, have a look at this, page 33 and beyond: PDF document. Best regards from V.

Re: Any expert for polynomial roots here? Message #3 Posted by fhub on 8 Nov 2011, 8:02 a.m., in response to message #2 by Valentin Albillo Quote: In case you haven't already, have a look at this, page 33 and beyond: PDF document. Oh my god, a 11MB PDF document! ;-) Thanks, I'll have a look at it, Franz

Re: Any expert for polynomial roots here? Message #4 Posted by Alexander Oestert on 8 Nov 2011, 9:13 a.m., in response to message #3 by fhub Quote: Oh my god, a 11MB PDF document! ;-) Don't tell me you're still using a 14.4 modem... ;o

Re: Any expert for polynomial roots here? Message #5 Posted by fhub on 8 Nov 2011, 9:59 a.m., in response to message #4 by Alexander Oestert Quote: Don't tell me you're still using a 14.4 modem... ;o No, I didn't mean the download time but my old (and very weak) Win98 notebook that I'm using here: with about 70MB free RAM it even uses the Windows swapfile when I scroll through a 11MB PDF document! :-( But I just saw that I have already exactly the same HP-Journal (found it when I made my HP-71 emulator package), but my file is only 1MB compared to the version above with 11MB (nevertheless it has the same content).

Re: Any expert for polynomial roots here? Message #6 Posted by Namir on 8 Nov 2011, 7:57 a.m., in response to message #1 by fhub Hello Franz, I use Bairstow's method when dealing with polynomials with real coefficients, because you should get all the roots--real or complex. HP seems to favor the Laguerre method. As far as speed is concerned, I never thought that Bairstow is slow. Maybe you need to double check you source of equations and your implementation. Namir

Re: Any expert for polynomial roots here? Message #7 Posted by fhub on 8 Nov 2011, 8:18 a.m., in response to message #6 by Namir Ho Namir, Quote: I use Bairstow's method when dealing with polynomials with real coefficients, because you should get all the roots--real or complex. HP seems to favor the Laguerre method. Well, if I don't misunderstand what I've read then also Laguerre should give all (real+complex) roots!? The only problem could be that with real coefficients you always have pairs of complex conjugate solutions, and Laguerre doesn't seem to work very well if there are solutions with equal or similar absolute values (which is of course true for conjugates). Quote: As far as speed is concerned, I never thought that Bairstow is slow. Maybe you need to double check you source of equations and your implementation. I guess you're talking about PC implementations!? Well, for a PC program the speed (of Bairstow) is of course no big deal - even a few 1000 iterations don't need much time - but on slow hardware like 20/30b the speed is indeed a very important factor. And this is the reason why I was searching for eventually better (=faster) methods than Bairstow (although less than 5 sec for a 10th-degree polynomial is acceptable IMO). I don't think that something in my implementation is wrong (or inefficient), both algorithms are so simple that it's almost impossible to do anything more complicated than necessary. Laguerre requires to do all calculations complex (which makes it a bit slower), but this should be compensated by what is stated in this WikiPedia article: "is only very rarely required to compute more than a few iterations to get high accuracy". Well, with Laguerre I get about 200-300 iterations for completely solving a 10-degree polynomial, compared to 60-70 with Bairstow. (?) Franz

Re: Any expert for polynomial roots here? Message #8 Posted by Artur-Brazil on 9 Nov 2011, 5:39 a.m., in response to message #7 by fhub Great algorithms for implementing in HP-15C !!!

Re: Any expert for polynomial roots here? Message #9 Posted by hugh steers on 9 Nov 2011, 8:48 p.m., in response to message #1 by fhub i can report that i use matrix methods for polynomial roots. essentially, i translate it into an eigenvalue problem. also i get complex results this way. results are usually stable and quite fast. im using a 14Mhz calculator for this. with this method, you dont have to extract pairs of roots at a time and reduce. so you dont have a forward error problem. however, Namir has done an in-depth study of polynomial root problems. you might want to try his examples. my method also has problems with multiple roots, where i get imprecise results. so my approach is not foolproof.

Re: Any expert for polynomial roots here? Message #10 Posted by fhub on 10 Nov 2011, 5:29 a.m., in response to message #9 by hugh steers Quote: i can report that i use matrix methods for polynomial roots. essentially, i translate it into an eigenvalue problem. also i get complex results this way. Yes, I know this method, but there are 2 problems (on the WP34s): First this would limit the polynomial solver to degree 10 (max. 10x10 matrices possible) and then the WP34s has no eigenvalues function, so _this_ would have to be written instead of a direct polyroots. And when I look at the method(s) for calculating eigenvalues (e.g. QR decompossition), I'm not sure if this would even fit into one flash region. Quote: my method also has problems with multiple roots, where i get imprecise results. so my approach is not foolproof. Well, multiple roots are of course no principle problem because they can always be removed by reducing the polynomial P to P/gcd(P,P'), but that would also need many additional program steps for getting the gcd(P,P'). After having read through all the infos about the Laguerre method, now I know what makes this polynomial solver on the HP-71 (which is based on the good old ZERPOL Fortran routine) so powerful: it's not the Laguerre method in the first, but the clever chosing and changing of initial values and/or iteration steps! This method calculates 4 different boundary values where one of them even needs the calculation of the roots of a slightly modified polynomial of the same degree (again with an iterative method: Newton-Raphson)! That means: just to get a good initial guess (or iteration step), you first have to solve almost the same polynomial with a different method, and then you can use one of these found roots (the minimum of 4 boundaries) for the Laguerre iteration! Well, that may be no problem for a PC (or a calculator with enough RAM), but it's definitely not suited for a programmable calculator. So in fact this Laguerre method isn't that good at all, it works only well if you have VERY good initial values in every iteration step. This is probably the explanation for my bad experiences with my Laguerre program on the WP34s compared to the Bairstow method. Franz Edited: 10 Nov 2011, 5:32 a.m.

Re: Any expert for polynomial roots here? Message #11 Posted by Valentin Albillo on 10 Nov 2011, 6:42 a.m., in response to message #10 by fhub Quote: Well, multiple roots are of course no principle problem because they can always be removed by reducing the polynomial P to P/gcd(P,P'), but that would also need many additional program steps for getting the gcd(P,P'). Nope. Computing the gcd of two polynomials is pretty straightforward and fast, so it can be done if deemed useful. Quote: This method calculates 4 different boundary values where one of them even needs the calculation of the roots of a slightly modified polynomial of the same degree (again with an iterative method: Newton-Raphson)! Nope. You misunderstood Fig. 4 in page 35 in the HP Journal issue I referred you to. It doesn't say "roots", plural, but "root", singular, i.e., you need just one root of the modified polynomial, which the Newton-Raphson method gets you in a flash using the very good initial estimate provided. So it isn't that hard at all. Quote: Well, that may be no problem for a PC (or a calculator with enough RAM), but it's definitely not suited for a programmable calculator. Nope. It is perfectly suitable, if correctly and efficiently implemented. Quote: So in fact this Laguerre method isn't that good at all, it works only well if you have VERY good initial values in every iteration step. This is probably the explanation for my bad experiences with my Laguerre program on the WP34s compared to the Bairstow method. Nope on both counts. The Laguerre method is very good, that's why it was selected as the basis for the PROOT keyword in the HP-71 Math ROM, and I think the explanation for your bad experiences lies somewhere else. Best regards from V.

Re: Any expert for polynomial roots here? Message #12 Posted by fhub on 10 Nov 2011, 7:09 a.m., in response to message #11 by Valentin Albillo Quote: Nope. Computing the gcd of two polynomials is pretty straightforward and fast, so it can be done if deemed useful. Yes, straightforward and fast, but you need 3 routines for P', for gcd and polynomial division, and that's not done in a few bytes. Quote: Nope. It is perfectly suitable, if correctly and efficiently implemented. Well, then why don't you make such a WP34s implementation? ;-) Quote: Nope on both counts. The Laguerre method is very good, that's why it was selected as the basis for the PROOT keyword in the HP-71 Math ROM, and I think the explanation for your bad experiences lies somewhere else. And what should this reason be? It can't be any error in my code because else I won't get correct results. The reason is probably that I just use the origin (0/0) as initial guess and then random guesses in the interval [-1,1] after every 25 unsuccessful iterations, but if the Laguerre method would really be as good as stated everywhere, than it would also work with _these_ guesses (and not require complicated calculations for 'good' initial values). Franz

Re: Any expert for polynomial roots here? Message #13 Posted by Valentin Albillo on 10 Nov 2011, 12:05 p.m., in response to message #12 by fhub Quote: Yes, straightforward and fast, but you need 3 routines for P', for gcd and polynomial division, and that's not done in a few bytes. Nope. It can be done in a few bytes if correctly and efficiently implemented. Quote: Well, then why don't you make such a WP34s implementation? ;-) Got much better things to do with my time, and besides "been there, done that, bought the T-shirt". I already implemented a full-fledged, efficient polynomial solver for the HP-41C back then in the golden days some 30 years ago so no need to repeat myself, see the relevant PPC issue if interested. Quote: And what should this reason be? It can't be any error in my code because else I won't get correct results. Nope, I haven't inspected your code but I'd hazard that the reason is your lack of relevant theoretical knowledge about the task you've set yourself up to. Quote: [...] but if the Laguerre method would really be as good as stated everywhere, than it would also work with _these_ guesses (and not require complicated calculations for 'good' initial values). Nope, they aren't complicated at all, it's your lack of understanding that makes them seem complicated to your eyes, thus resulting in you unfairly critizicing something which you simply fail to understand. Quote: I've read again this part, but I'm still not sure who of us is right!? I'll give you a hint: you're wrong. Quote: [...] so IMO to be sure to really get the minimum you would have to calculate all roots. Nope, you don't need all roots at all, you only need a very specific one. It's just your lack of relevant theoretical understanding getting in the way, may I suggest you get and study any of the many good books on the subject, "Higher Algebra" by Kurosh is one of the best. Quote: So as I already stated: the 'power' of this PROOT routine is not the Laguerre method, but the big effort that is put into finding a good initial guess. But with such a good initial guess almost every other polynomial rootfinding algorithm would work just as well. Nope, you're wrong in all counts but I won't push the matter further, I've already suggested that you learn the theory of what you're trying to accomplish lest you'll waste tons of effort to get a mediocre implementation at best. Best regards from V.

Re: Any expert for polynomial roots here? Message #14 Posted by fhub on 10 Nov 2011, 12:30 p.m., in response to message #13 by Valentin Albillo It seems "nope" is the only word in your quite limited vocabulary! And I'm really tired of discussing with such an arrogant and unfriendly guy. Unfortunately you're not the only one of this kind here in this forum, and so I'm now finally leaving this place full of all-knowing 'gods'. Bye, Franz.

Re: Any expert for polynomial roots here? Message #15 Posted by Walter B on 11 Nov 2011, 2:56 a.m., in response to message #14 by fhub Quote: ... so I'm now finally leaving this place full of all-knowing 'gods'. 17 and counting [;-)

Re: Any expert for polynomial roots here? Message #16 Posted by Bart (UK) on 11 Nov 2011, 5:02 a.m., in response to message #14 by fhub Quote: this place full of all-knowing 'gods'. Of which you are one as well. I probably am too at times :-)

Re: Any expert for polynomial roots here? Message #17 Posted by Ángel Martin on 10 Nov 2011, 3:26 p.m., in response to message #13 by Valentin Albillo Here´s Valentin´s program from PPC TN - slightly changed for the 41Z consumption (very obvious places and easy to undo) 1 LBL "ZPROOT" 87 E-3 2 SIZE? 88 ST+ 01 3 "DEGREE=?" 89 RCL 03 4 PROMPT 90 STO IND 05 5 STO Z 91 RCL 04 6 ST+X 92 STO IND 06 7 11 93 DSE 00 8 + 94 GTO 06 9 X>Y? 95 TONE 5 10 PSIZE 96 RCL 01 11 RCL Z 97 INT 12 STO 00 98 E1 13 STO 03 99 - 14 9,008 100 E3 15 + 101 / 16 STO 01 102 ST- 05 17 STO 05 103 FIX 3 18 X<>Y 104 SF 21 19 E 105 LBL 10 20 - 106 ISG 00 21 STO 02 107 NOP 22 STO 06 108 RCL IND 06 23 FIX 0 109 RCL IND 05 24 CF 29 110 ZAVIEW 25 LBL 05 111 DSE 06 26 "IM^RE(" 112 DSE 05 27 ARCL 03 113 GTO 10 28 "@)=?" 114 CF 21 29 PROMPT 115 SF 29 30 STO IND 05 116 RTN 31 X<>Y 117 LBL 11 32 STO IND 06 118 RCL 01 33 DSE 03 119 STO 05 34 X<>Y 120 RCL 02 35 DSE 06 121 STO 06 36 DSE 05 122 FC? 01 37 GTO 05 123 GTO 13 38 RCL 03 124 E-3 39 LBL 06 125 ST+ 05 40 "SOLVING..." 126 LBL 13 41 AVIEW 127 RCL IND 06 42 SF 25 128 RCL IND 05 43 SF 99 129 FC? 01 44 CF 00 130 GTO 02 45 CHS 131 RCL 08 46 STO 04 132 ST* Z 47 FIX 2 133 * 48 RND 134 DSE 08 49 FIX 6 135 GTO 02 50 X#0? 136 RTN 51 GTO 01 137 LBL 00 52 SIGN 138 ZENTER^ 53 STO 04 139 RCL 04 54 LBL 01 140 RCL 03 55 RCL 00 141 Z* 56 STO 08 142 RCL IND 05 57 SF 01 143 FS? 01 58 XEQ 11 144 RCL 08 59 R-P 145 FS? 01 60 1/X 146 * 61 STO 07 147 + 62 X<>Y 148 FS? 00 63 CHS 149 STO IND 05 64 STO 08 150 X<>Y 65 CF 01 151 RCL IND 06 66 XEQ 11 152 FS? 01 67 ZENTER^ 153 RCL 08 68 RCL 08 154 FS? 01 69 RCL 07 155 * 70 P-R 156 + 71 Z* 157 FS? 00 72 ST- 03 158 STO IND 06 73 X<>Y 159 X<>Y 74 ST- 04 160 FS? 01 75 ZRND 161 DSE 08 76 Z#0? 162 LBL 02 77 GTO 01 163 DSE 06 78 FIX 0 164 DSE 05 79 "FOUND ROOT#" 165 GTO 00 80 ARCL 00 166 END 81 AVIEW 82 SF 00 83 XEQ 11 84 E 85 ST+ 05 86 ST+ 06

Re: Any expert for polynomial roots here? Message #18 Posted by Ángel Martin on 11 Nov 2011, 2:12 a.m., in response to message #17 by Ángel Martin and here are the roots found by the program for Valentin´s second example: x^20+10 x^19+55 x^18+210 x^17+615 x^16+1452 x^15+2850 x^14+ 4740 x^13+6765 x^12+8350 x^11+8953 x^10+8350 x^9+6765 x^8+ 4740 x^7+2850 x^6+1452 x^5+615 x^4+210 x^3+55 x^2+10 x+1 = 0 1,473-J7,264 1,473+J7,264 0,476-J1,533 0,476+J1,533 -0,769-J0,180 -0,769+J0,180 -0,593-J0,367 -0,593+J0,367 0,138-J1,013 0,138+J1,013 -0,028-J0,805 -0,028+J0,805 -0,146-J0,684 -0,146+J0,684 -0,247-J0,599 -0,247+J0,599 -0,346-J0,528 -0,346+J0,528 -0,456-J0,458 -0,456+J0,458 execution time on V41: about 8 seconds (TURBO mode of course). It´s included in the 41Z module, so just load the image MOD file and turn it loose. Write or wrong? Easy enough to check ... Edited: 11 Nov 2011, 2:13 a.m.

Re: Any expert for polynomial roots here? Message #19 Posted by Gerson W. Barbosa on 11 Nov 2011, 10:20 a.m., in response to message #18 by Ángel Martin According to Wolfram|Alpha there are the following coincident roots: x = -1/2 - sqrt(3)/2*i x = -1/2 + sqrt(3)/2*i Cheers, Gerson.

Re: Any expert for polynomial roots here? Message #20 Posted by fhub on 10 Nov 2011, 10:26 a.m., in response to message #11 by Valentin Albillo Quote: Nope. You misunderstood Fig. 4 in page 35 in the HP Journal issue I referred you to. It doesn't say "roots", plural, but "root", singular, i.e., you need just one root of the modified polynomial, which the Newton-Raphson method gets you in a flash using the very good initial estimate provided. So it isn't that hard at all. I've read again this part, but I'm still not sure who of us is right!? This Cauchy lower bound RA is defined there as "minimum positive zero" of this modified polynomial, so IMO to be sure to really get the minimum you would have to calculate all roots. Ok, with the given (small) initial guess you'll rather get the minimum root first, but it's not absolutely sure. And furthermore I wonder why I would need this Cauchy bound as a 4th guess, when the 3 previous guesses are already so good!? So as I already stated: the 'power' of this PROOT routine is not the Laguerre method, but the big effort that is put into finding a good initial guess. But with such a good initial guess almost every other polynomial rootfinding algorithm would work just as well. Franz. Edited: 10 Nov 2011, 10:27 a.m.

Re: Any expert for polynomial roots here? Message #21 Posted by Peter Murphy (Livermore) on 10 Nov 2011, 1:52 p.m., in response to message #1 by fhub I am shocked.

[ Return to Index | Top of Index ]

Go back to the main exhibit hall