# HP Forums

Full Version: Any expert for polynomial roots here?
You're currently viewing a stripped down version of our content. View the full version with proper formatting.

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

In case you haven't already, have a look at this, page 33 and beyond:

Best regards from V.

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

Quote:

In case you haven't already, have a look at this, page 33 and beyond:

Oh my god, a 11MB PDF document! ;-)

Thanks, I'll have a look at it,

Franz

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

Quote:

Oh my god, a 11MB PDF document! ;-)

Don't tell me you're still using a 14.4 modem... ;o

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).

Great algorithms for implementing in HP-15C !!!

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.

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.

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.

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

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.

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.

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.

I am shocked.

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
```

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.

Quote:
... so I'm now finally leaving this place full of all-knowing 'gods'.

17 and counting [;-)

Quote:
this place full of all-knowing 'gods'.

Of which you are one as well. I probably am too at times :-)

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.