WP34S: Valentin Albillo's Polynomial Root Finder - Printable Version +- HP Forums (https://archived.hpcalc.org/museumforum) +-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html) +--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html) +--- Thread: WP34S: Valentin Albillo's Polynomial Root Finder (/thread-205064.html) WP34S: Valentin Albillo's Polynomial Root Finder - Miguel Toro - 11-17-2011 Hi, Here is a version for the WP 34s of Albillo's original program for the HP 41. I would like to thank Ángel Martin for helping me with lot of information, documentation and his own version for the 41z. Porting the program was very straightforward. The only real difference I found between commands was with the SIGN function that behaves like it does in the 3Xs family rather than it does in the 41 and 42 calculators, so I have to replace it in my code. The program has been tested with version 2.2 rev. 1782. Just copy it in a file and assemble it to try in the emulator and the calculator Testing the examples from the article on the calculator, I obtained the same results, with the following time: Example 1) 5 seconds – original on HP 41: 8 minutes : (2+8i)z^6 + (3+0i)z^5 + (-1+2i)z^4 + (0+2i)z^3 + (-3-3i)z^2 + (1+2i)z + (-2+3i) =0 Example 2) 3,5 seconds – original on HP 41: 5 minutes : 5x^6 – 4x^5 – 3x^4 + 8x^3 + 8x^2 – 2x + 7 = 0 Note: You can delete the STOP in line 34. It is there just for testing. Enjoy! Miguel ```001 LBL'PRT' 002 SSIZE4 003 FIX 00 004 CL[alpha] 005 [alpha]'N?' 006 PROMPT 007 STO 00 008 STO 03 009 STO+ X 010 9 011 . 012 0 013 0 014 8 015 0 016 2 017 + 018 STO 01 019 STO 05 020 CLx 021 LBL 05 022 CL[alpha] 023 [alpha]'IM' 024 [alpha] [^] 025 [alpha]'RE[space]' 026 [alpha]IP 03 027 [alpha] ? 028 PROMPT 029 [cmplx]STO[->]05 030 DEC 03 031 DSE 05 032 GTO 05 033 RCL 03 034 STOP 035 LBL 06 036 CF A 037 +/- 038 STO 04 039 FIX 02 040 ROUND 041 FIX 06 042 x[!=]0? 043 GTO 01 044 EEX 045 STO 04 046 LBL 01 047 RCL 00 048 STO 08 049 SF 01 050 XEQ 11 051 [->]POL 052 1/x 053 STO 07 054 x[<->] Y 055 +/- 056 STO 08 057 CF 01 058 XEQ 11 059 [cmplx]ENTER 060 [cmplx]RCL 07 061 [->]REC 062 [cmplx][times] 063 [cmplx]STO- 03 064 [cmplx]ROUND 065 [cmplx]x[!=]0? 066 GTO 01 067 SF A 068 XEQ 11 069 2 070 STO+ 05 071 EEX 072 3 073 +/- 074 [times] 075 STO+ 01 076 RCL 04 077 RCL 03 078 [cmplx]STO[->]05 079 x[<->] Y 080 DSE 00 081 GTO 06 082 RCL 01 083 IP 084 1 085 0 086 - 087 EEX 088 3 089 / 090 STO- 05 091 FIX 07 092 LBL 10 093 INC 00 094 CL[alpha] 095 [alpha]'RT[space]' 096 [alpha]IP 00 097 [alpha]VIEW 098 PSE 15 099 CL[alpha] 100 [alpha]'RE' 101 [cmplx]RCL[->]05 102 PROMPT 103 CL[alpha] 104 [alpha]'IM' 105 x[<->] Y 106 PROMPT 107 DSE 05 108 GTO 10 109 CLx 110 CF A 111 RTN 112 LBL 11 113 RCL 01 114 STO 05 115 FC? 01 116 GTO 13 117 EEX 118 3 119 +/- 120 ENTER[^] 121 2 122 [times] 123 STO+ 05 124 LBL 13 125 [cmplx]RCL[->]05 126 FC? 01 127 GTO 02 128 RCL 08 129 STO[times] Z 130 [times] 131 DSE 08 132 GTO 02 133 RTN 134 LBL 00 135 [cmplx]ENTER 136 [cmplx]RCL 03 137 [cmplx][times] 138 RCL[->]05 139 FS? 01 140 RCL 08 141 FS? 01 142 [times] 143 + 144 FS? A 145 STO[->]05 146 x[<->] Y 147 INC 05 148 RCL[->]05 149 FS? 01 150 RCL 08 151 FS? 01 152 [times] 153 + 154 FS? A 155 STO[->]05 156 DEC 05 157 x[<->] Y 158 FS? 01 159 DSE 08 160 LBL 02 161 DSE 05 162 GTO 00 163 RTN ``` Edited: 17 Nov 2011, 1:06 p.m. Re: WP34S: Valentin Albillo's Polynomial Root Finder - Ángel Martin - 11-18-2011 Bravo Miguel! This contributes to create a 34S library :-) The next should be the Bessel functions - Complex of course - sure the 41Z programs can also help! Looking at the 34S code I can recognize a couple of tricks I used to adapt the program to the 41Z - nice to see the design similarities, which on the other hand are pretty logical choices. I'm thinking about Valentin's examples in the other thread - would be good to test the port with those and see how it fares, sort of "turning the creature against the creator" :-) the first one was: ```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 ``` which Wolfram Alpha reports (Thanks Gerson!) having coincident roots at: ```x = -1/2 - sqrt(3)/2*i x = -1/2 + sqrt(3)/2*i ``` Not quite the same as what this program did when I tried it on V41/ w/41Z - I may have made some data entry error so I'll repeat it this week-end. Best, 'AM Re: WP34S: Valentin Albillo's Polynomial Root Finder - Paul Dale - 11-18-2011 Quote:Looking at the 34S code I can recognize a couple of tricks I used to adapt the program to the 41Z - nice to see the design similarities, which on the other hand are pretty logical choices. What in particular? - Pauli Re: WP34S: Valentin Albillo's Polynomial Root Finder - Ángel Martin - 11-18-2011 The way to push the complex argument into the complex stack with independence from the real one, and also checking on the complex round as opposed to on each real and imaginary parts - as it was done in the original program. Trivialities when seen in retrospective but it was different when you also had to write the MCODE underneath all the functions to play well in that scheme. ÁM Re: WP34S: Valentin Albillo's Polynomial Root Finder - Paul Dale - 11-18-2011 Quote:The way to push the complex argument into the complex stack with independence from the real one ... We don't separate the stack into a complex and real sections -- everything lives on the single stack. X + iY, Z + iT, A + iB and C + iD. The latter two being optional. I did add complex STO and RCL commands to make life easier. - Pauli Edited: 18 Nov 2011, 2:36 a.m. Re: WP34S: Valentin Albillo's Polynomial Root Finder - Dieter - 11-18-2011 Looking at the code, I would like to add a few comments: I would prefer the input routine to prompt for the coefficients with a zero in X. This way real coefficents can be entered much easier, without the need of typing [0] [ENTER] first. Is there a reason why line 117-122ff do not simply contain a straightforward constant ". 0 0 2" or "2 EEX +/- 3" ?-) Also, line 071ff. muliplies by 1E-3. Why not simply divide by 1E3 ? Unlike the 41-series, the 34s can do RCL-Arithmetics. So step 150 - 154 can be replaced by a simple FS?01  RCLx 08. Is the complex ENTER in line 059 and 135 really required? The program switches between four different display formats (FIX 00, 02, 06 and 07). Why is this neccessary? I see some ROUND statements here and there. What exactly do they do? The ROUND in step 064 has no preceding display setting so it's unclear how many digits are rounded. BTW, there also is a "X~~0" test that does not require rounding at all. I assume the algorithm iterates as long as a certain precision is not yet achieved, and the test in step 064 checks this. In FIX 06 mode this means the absolute error is less than 5E-7. This is fine for small and moderate values, but then think of a root near 1E-6. Or another one near 1E9. So there might be a better way: a root is found if the latest and the previous value agree in 12 significant digits. Which should be fine since the 34s works with 16 digits. The EEX in step 044 probably is a leftover from an earlier 41-implemenation where it was a few milliseconds faster than a usual "1" (one). I don't think this is required on the 34s. ;-) You see there is some nitpicking in these comments. Please consider this a compliment since there obviously are no major bugs. ;-) Dieter Re: WP34S: Valentin Albillo's Polynomial Root Finder - Miguel Toro - 11-18-2011 Hi Dieter, Excellent and pertinent remarks! There is a lot of heritage from the original code that is not required with the WP 34s. In this version I just made the more obvious changes and improvements, thanks to the excellent complex instruction set of the 34S, just to test how easily was to port code from the 41. I will modify the program based on your comments. Thanks and regards, Miguel Re: WP34S: Valentin Albillo's Polynomial Root Finder - Dieter - 11-18-2011 Here's a first attempt of a modified version for the 34s. I'm sure there's room for improvements, maybe there are even some errors. The algorithm is the same, basically the following details have been changed: The various rounding steps have been removed. Instead, a tolerance is stored in register A, here with a value of 1E-6. From this, three individual "deltas" are derived: The original test if abs(R04) < 0,005 was replaced by R04^2 < 1E-6. In the output routine, roots with imaginary parts less than 1E-12 are neglected, i.e. the root is considered real. And finally there is a more sophisticated convergence check: First, the real parts are compared. If they agree in all displayed digits (i.e. all 12, since we're in ALL mode), the imaginary parts are compared as well. If they do not agree the program finally tests whether both imaginary parts are smaller than approx. 1E-12. In this case the root is virtually real and no further iteration is required. The input routine now prompts for the coefficients with a cleared stack. So real values can be entered without pressing any additional keys. The output routine has been rewritten. Every root is now labelled as "Real" or "Compl". In the latter case the result is x+iy and both parts can be viewed with x<>y, as usual. Finally, various minor improvements have been applied. ```001 LBL'PRT' 002 SSIZE4 003 CLSTK 004 CF 01 005 CF A 006 CL[alpha] 007 [alpha]'N=?' 008 PROMPT 009 x>0? 010 FP? 011 BACK 03 012 STO 00 013 STO 03 014 STO+ X 015 9 016 . 017 0 018 0 019 8 020 0 021 2 022 + 023 STO 01 024 STO 05 025 LBL 05 026 CL[alpha] 027 [alpha]'Im' 028 [alpha] [^] 029 [alpha]'Re[space]' 030 [alpha]IP 03 031 [alpha] ? 032 CLSTK 033 PROMPT 034 [cmplx]STO[->]05 035 DEC 03 036 DSE 05 037 GTO 05 038 EEX 039 +/- 040 6 041 STO A // tolerance = 1E-6 042 ALL 03 043 RCL 03 044 LBL 06 045 CF A 046 +/- 047 STO 04 048 x[^2] 049 x[>=]? A // i.e. |R04| >= 1E-3? 050 GTO 01 051 1 052 STO 04 053 LBL 01 054 RCL 00 055 STO 08 056 SF 01 057 XEQ 11 058 [->]POL 059 1/x 060 STO 07 061 x[<->] Y 062 +/- 063 STO 08 064 CF 01 065 XEQ 11 066 [cmplx]RCL 07 067 [->]REC 068 [cmplx][times] 069 [cmplx]+/- 070 [cmplx]RCL+ 03 071 [cmplx]ENTER 072 [cmplx]x[<->] 03 073 x[approx]? Z 074 SKIP 01 075 GTO 01 076 x[<->] T 077 x[approx]? Y 078 GTO 04 079 [cmplx]ABS 080 [sqrt] 081 x>? A // both imaginary parts < 1E-12? 082 GTO 01 083 LBL 04 084 SF A 085 XEQ 11 086 2 087 STO+ 05 088 EEX 089 3 090 / 091 STO+ 01 092 RCL 04 093 RCL 03 094 [cmplx]STO[->]05 095 x[<->] Y 096 DSE 00 097 GTO 06 098 LBL'RES' 099 SF A 100 RCL 01 101 FP 102 9 103 + 104 STO 05 105 CLx 106 STO 00 107 LBL 10 108 INC 00 109 CL[alpha] 110 [alpha]'Rea' 111 [alpha]'l[space]' 112 [cmplx]RCL[->]05 113 x[<->] Y 114 ABS 115 x]05 122 PROMPT 123 ISG 05 124 GTO 10 125 CLSTK 126 CF A 127 RTN 128 LBL 11 129 . 130 0 131 0 132 2 133 RCL 01 134 FS? 01 135 + 136 STO 05 137 [cmplx]RCL[->]05 138 FC? 01 139 GTO 02 140 RCL 08 141 STO[times] Z 142 [times] 143 DSE 08 144 GTO 02 145 RTN 146 LBL 00 147 [cmplx]RCL[times] 03 148 RCL[->]05 149 FS? 01 150 RCL[times] 08 151 + 152 FS? A 153 STO[->]05 154 x[<->] Y 155 INC 05 156 RCL[->]05 157 FS? 01 158 RCL[times] 08 159 + 160 FS? A 161 STO[->]05 162 DEC 05 163 x[<->] Y 164 FS? 01 165 DSE 08 166 LBL 02 167 DSE 05 168 GTO 00 169 RTN ``` Here's an example: x^4 - 4 x^3 - 3 x^2 + x + 1 ```XEQ"PRT" N=? 4 [R/S] Im^Re 4? 1 [R/S] Im^Re 3? -4 [R/S] Im^Re 2? -3 [R/S] Im^Re 1? 1 [R/S] Im^Re 0? 1 [R/S] Compl 1 = -0,57898166378 [x<>y] 0,22687379984 [R/S] Compl 2 = -0,57898166378 [x<>y] -0,22687379984 [R/S] Real 3 = 0,56277097739 [R/S] Real 4 = 4,59519235018 [R/S] 0 ``` Once calculated, the results can be recalled by XEQ"RES". By the way: Walter, I miss a directly accessible "=" on the alpha keyboard. As far as I can see, common prompts like the "N=?" in line 007 cannot be entered from the keyboard. However, there are several unassigned keys on the alpha keyboard: Just as h [4] enters a question mark, h [5] could be used for the missing "=". Dieter Re: WP34S: Valentin Albillo's Polynomial Root Finder - Walter B - 11-18-2011 Dieter, thanks for pointing this out. I'd put '=' on h-shifted ENTER in alpha mode - the location CONST lives. Re: WP34S: Valentin Albillo's Polynomial Root Finder - Marcus von Cube, Germany - 11-18-2011 '=' is in the test catalog of course. In multi alpha mode I've put it on the same key (no catalogs available then) in the most recent build. h ENTER^ produces an up arrow. Re: WP34S: Valentin Albillo's Polynomial Root Finder - Walter B - 11-19-2011 As often before, this is not finally discussed nor agreed on yet. Some procedures seem quite difficult to teach ;-) And a catalogue access is not a direct access 8-) Edit: Meanwhile, we agreed on putting '=' on h-shifted R/S in alpha mode, the key it lives nearly everywhere ;-) Edited: 19 Nov 2011, 5:15 p.m. Re: WP34S: Valentin Albillo's Polynomial Root Finder - fhub - 11-20-2011 Quote: Here's a first attempt of a modified version for the 34s. I'm sure there's room for improvements, maybe there are even some errors. Hi Dieter, it seems you have indeed changed it a bit too much - your version doesn't even solve the simple polynomial x^10+1=0. No error messages, no wrong solution(s), just runs forever ... Franz Re: WP34S: Valentin Albillo's Polynomial Root Finder - fhub - 11-20-2011 Quote: I'm thinking about Valentin's examples in the other thread - would be good to test the port with those and see how it fares, sort of "turning the creature against the creator" :-) the first one was: ```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 ``` which Wolfram Alpha reports (Thanks Gerson!) having coincident roots at: ```x = -1/2 - sqrt(3)/2*i x = -1/2 + sqrt(3)/2*i ``` Not quite the same as what this program did when I tried it on V41/ w/41Z - I may have made some data entry error so I'll repeat it this week-end. Best, 'AM Hi Angel, either you have mistyped any coefficient(s), or Miguel's translation is not correct, or Valentin's program does indeed not solve this polynomial at all!? I've entered this polynomial now 3 or 4 times and waited about 2 minutes (on the PC emulator of WP34s), but no solution - it just didn't stop running. So I guess I'll stay with my own new polynomial solver (with Laguerre method), which solves all those polynomials without any problems (even the 2 'hard' ones that Valentin has posted), and is in fact much faster and gives more accurate solutions. :-) Franz Re: WP34S: Valentin Albillo's Polynomial Root Finder - Dieter - 11-20-2011 Yes, the more sophisticated convergence test was a bit too ...sophisticated. ;-) It does not work if a root is very close to zero and decides that two consecutive approximations do not agree even if their real parts are 1E-24 and 1E-25. In your example the problem arises at the fifth and sixth root (i and -i). Their real part is exactly zero while the program actually has determined a value near 1E-16. Here's an improved version. It checks whether the real and complex parts either agree in all displayed digits or if they both are virtually zero. ```... ... 058 [->]POL 059 1/x 060 [cmplx]CONJ // shorter than 061 [cmplx]STO 07 // previous code 062 STO 07 063 CF 01 064 XEQ 11 065 [cmplx]RCL 07 066 [->]REC 067 [cmplx][times] 068 [cmplx]+/- 069 [cmplx]RCL+ 03 070 [cmplx]ENTER 071 [cmplx]x[<->] 03 072 x[<->] Y 073 R[v] // x=re_old y=re_new z=im_new z=im_old 074 x[approx]? Y // do both real parts agree? 075 SKIP 04 // then continue checking 076 [cmplx]ABS // both real parts do not agree 077 [sqrt] // so check whether they are both near zero 078 x>? A // are both real parts < 1E-12? 079 GTO 01 // if not, do next iteration 080 R[v] // real parts agree or both are near zero 081 R[v] // now check imaginary parts 082 x[approx]? Y // do both imaginary parts agree as well? 083 SKIP 04 // then we're done 084 [cmplx]ABS // at this point both real parts agree 085 [sqrt] // but the imaginary parts don't 086 x>? A // now check if they are at least < 1E-12 087 GTO 01 // if not, do another iteration 088 SF A // otherwise we're done, continue here 089 XEQ 11 ... ... ``` Dieter Edited: 20 Nov 2011, 9:12 a.m. Re: WP34S: Valentin Albillo's Polynomial Root Finder - fhub - 11-20-2011 Quote: Yes, the more sophisticated convergence test was a bit too ...sophisticated. ;-) Have you tried Valentin's 20-degree polynomial which Angel has re-posted above? Let's see if your new version is now sophisticated enough for this one ... ;-) Or you could even try the 30-degree polynomial Valentin has posted as test for my first (Bairstow) solver. It has the only solution x=-1 (but with multiplicity 30), and my new Laguerre solver gives the exact solution 30-times! :-) Franz Re: WP34S: Valentin Albillo's Polynomial Root Finder - Ángel Martin - 11-20-2011 My questions exactly :-) Re: WP34S: Valentin Albillo's Polynomial Root Finder - fhub - 11-20-2011 Quote: My questions exactly :-) Well, I've tried at least the 20-degree poly (with Miguel's translated program),but got no solution at all. So I didn't try the 30-degree one. Re: WP34S: Valentin Albillo's Polynomial Root Finder - Dieter - 11-20-2011 Franz, Quote: Let's see if your new version is now sophisticated enough for this one ... ;-) there obviously is a misunderstanding here. This is not "my" solution of a polynomial root solver. Honestly, I do not plan on writing something like this - my knowledge on this topic is very limited. What you see is simply Miguel's original PRT version with some minor improvements (input/output, convergence test, special 34s features). The whole algorithm is exactly the same as before. So if the original version does not handle the 20-degree-polynomial, the modified version will not do so either. So it's about time for you to show your latest Laguerre solver. ;-) Dieter Edited: 20 Nov 2011, 2:08 p.m. Re: WP34S: Valentin Albillo's Polynomial Root Finder - fhub - 11-20-2011 Quote: So it's about time for you to show your latest Laguerre solver. ;-) No Dieter, not after these attacks a week ago! Since this time my programs are only private. Franz Re: WP34S: Valentin Albillo's Polynomial Root Finder - Andrés C. Rodríguez (Argentina) - 11-21-2011 Franz, Of course, it's your choice to share your knowledge with the community or not; but it's puzzling that such a decision is made upon the result of a single event or exchange with one or two persons. Many contributors, even very bright ones, received some unexpected, unfair or out-of-style feedback from time to time; I recall some cases and I'm sure there are more. And sometimes the language differences, personal styles, local/country jargon, sense of irony, and translations add up to make a supossedly friendly message look more aggressive than intended. Sometimes any of us may be on the receiving side of such messages and sometimes, even by accident, may be at the originator side. But the decision to share or not is up to each one of us. Many people come here just to learn from others' creations and discoveries, and an implicit thank you is available for all contributors each time they post interesting and useful material. Please see that I'm not trying to analyze, explain, or otherwise support or condemn the opinions and postings of any other member of this community. Best regards. (I apologize in advance if my imperfect command of the English language prevents this message to convey its friendly and positive intentions) Re: WP34S: Valentin Albillo's Polynomial Root Finder - Walter B - 11-21-2011 Just crossed my mind whether we need a general disclaimer all of us writing in English as a foreign language should append to our postings in this forum. This won't, however, cover the numerous spelling or even wording errors by authors with "English" as their mother tongue observed here. Anyway, since it won't take care of all cases I remember having seen here, forget it ;-) Re: WP34S: Valentin Albillo's Polynomial Root Finder - Ángel Martin - 11-21-2011 what about this for a disclaimer: "get a grip" :-) Re: WP34S: Valentin Albillo's Polynomial Root Finder - Miguel Toro - 11-22-2011 Quote: the first one was: 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 which Wolfram Alpha reports (Thanks Gerson!) having coincident roots at: x = -1/2 - sqrt(3)/2*i x = -1/2 + sqrt(3)/2*i I tried this example with the original 41 program entered in the V41 emulator and it does not converge either. I did it several times in case I made any mistake entering the data. So I think that my translation is as faithful as it can be but one never knows... With both versions in FIX 6, results are really similar; and putting simply FIX 11 in line 041 or the wp34s version, accuracy is very much improved for many other cases I tested, but this one in particular, seems to be really tough :-). As disclaimed by Valentin in his article : "convergence is not guaranteed". What I am going to do: 1.- let the program as simple and near the original as possible, 2.- include some of the improvements suggested by Dieter, 3.- add an iteration limit so the program may stop with a “no zero found” message, 4.- put it in the Article section, so there is at least a solver available until some improved version is made by one of the many knowledgeable peoples of this community (someone please give a WP 34S to Valentin! I am sure he will enjoy). Thanks and regards, Miguel Edited: 22 Nov 2011, 10:35 a.m. Re: WP34S: Valentin Albillo's Polynomial Root Finder - Peter Murphy (Livermore) - 11-22-2011 Walter, you are correct: there are so many Forum participants who write excellent English despite its not being their mother tongue that it is easy for native speakers of English to assume that everyone here can communicate with perfect clarity. You are also correct that some cases seen here do not result from imperfect English: perfect clarity is sometimes the problem. My suggestion: if we treat our communications here as if they were delivered face to face, we might be less likely to speak in ways that give offense. And in cases of seeming offense, we might do well to let a matter rest overnight before responding, if we respond at all. In cases of offense, sincere and humble apology does wonders. Finally, remembering our purpose in participating in the Forum, namely to increase the enjoyment and utility we and others get from HP calculators, might do a lot to reduce the number of cases arising. Re: WP34S: Valentin Albillo's Polynomial Root Finder - Werner - 11-23-2011 I'd like to see Valentin's original program. I searched through PPC Journals in vain for it. BTW since his 20-degree polynomial is simply ` (x^2 + x + 1)^10 ` it's no wonder all polynomial solvers choke on it - well at least those that do not check for multiple roots up front (which can only really be done when you have exact coefficients, as is the case here). For a root with multiplicity n, (n-1)/n of the digits of the result will be inaccurate. Cheers, Werner Re: WP34S: Valentin Albillo's Polynomial Root Finder - Ángel Martin - 11-23-2011 Hi Werner, the reference is: PPC TECHNICAL NOTES (the Australian Chapter if I´m not mistaken) V1 N3 October 1980 Hope this helps, I can mail it to you if you still can't locate it - let me know. Best, 'AM Re: WP34S: Valentin Albillo's Polynomial Root Finder - Werner - 11-23-2011 Can't locate it.. it's not on TOS, and not on the Museum's DVD.. Werner Re: WP34S: Valentin Albillo's Polynomial Root Finder - Miguel Toro - 11-23-2011 Hi, You can find the program here in the article section. There is now an iteration limit in case there is no zero found. Regards, Miguel Edited: 23 Nov 2011, 7:40 p.m. Re: WP34S: Valentin Albillo's Polynomial Root Finder - Werner - 11-25-2011 Thanks for the file, Angel. I am, however, slightly disappointed as it is a simple Newton-Rhapson scheme, not an implementation of the Laguerre method, which I had expected. Cheers, Werner