Posts: 1,216
Threads: 75
Joined: Jun 2011
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 35 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
Posts: 1,755
Threads: 112
Joined: Jan 2005
In case you haven't already, have a look at this, page 33 and beyond:
PDF document.
Best regards from V.
Posts: 2,247
Threads: 200
Joined: Jun 2005
Hello Franz,
I use Bairstow's method when dealing with polynomials with real coefficients, because you should get all the rootsreal 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
Posts: 1,216
Threads: 75
Joined: Jun 2011
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
Posts: 1,216
Threads: 75
Joined: Jun 2011
Ho Namir,
Quote:
I use Bairstow's method when dealing with polynomials with real coefficients, because you should get all the rootsreal 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 10thdegree 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 200300 iterations for completely solving a 10degree polynomial, compared to 6070 with Bairstow. (?)
Franz
Posts: 429
Threads: 31
Joined: Jul 2011
Quote:
Oh my god, a 11MB PDF document! ;)
Don't tell me you're still using a 14.4 modem... ;o
Posts: 1,216
Threads: 75
Joined: Jun 2011
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 HPJournal (found it when I made my HP71 emulator package), but my file is only 1MB compared to the version above with 11MB (nevertheless it has the same content).
Posts: 127
Threads: 21
Joined: Jul 2006
Great algorithms for implementing in HP15C !!!
Posts: 536
Threads: 56
Joined: Jul 2005
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 indepth 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.
Posts: 1,216
Threads: 75
Joined: Jun 2011
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 HP71 (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: NewtonRaphson)!
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.
Posts: 1,755
Threads: 112
Joined: Jan 2005
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: NewtonRaphson)!
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 NewtonRaphson 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 HP71 Math ROM, and I think the explanation for your bad experiences lies somewhere else.
Best regards from V.
Posts: 1,216
Threads: 75
Joined: Jun 2011
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 HP71 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
Posts: 1,216
Threads: 75
Joined: Jun 2011
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 NewtonRaphson 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.
Posts: 1,755
Threads: 112
Joined: Jan 2005
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 Tshirt". I already implemented a fullfledged, efficient polynomial solver for the HP41C 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.
Posts: 1,216
Threads: 75
Joined: Jun 2011
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 allknowing 'gods'.
Bye,
Franz.
Posts: 167
Threads: 33
Joined: Jul 2011
Posts: 1,253
Threads: 117
Joined: Nov 2005
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 E3
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 E3
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 RP 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 PR 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
Posts: 1,253
Threads: 117
Joined: Nov 2005
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,473J7,264
1,473+J7,264
0,476J1,533
0,476+J1,533
0,769J0,180
0,769+J0,180
0,593J0,367
0,593+J0,367
0,138J1,013
0,138+J1,013
0,028J0,805
0,028+J0,805
0,146J0,684
0,146+J0,684
0,247J0,599
0,247+J0,599
0,346J0,528
0,346+J0,528
0,456J0,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.
Posts: 4,587
Threads: 105
Joined: Jul 2005
Quote:
... so I'm now finally leaving this place full of allknowing 'gods'.
17 and counting [;)
Posts: 850
Threads: 10
Joined: Mar 2009
Quote:
this place full of allknowing 'gods'.
Of which you are one as well. I probably am too at times :)
Posts: 2,761
Threads: 100
Joined: Jul 2005
According to WolframAlpha there are the following coincident roots:
x = 1/2  sqrt(3)/2*i
x = 1/2 + sqrt(3)/2*i
Cheers,
Gerson.
