The Museum of HP Calculators

HP Forum Archive 20

[ Return to Index | Top of Index ]

Polynomial roots again ...
Message #1 Posted by fhub on 23 Dec 2011, 4:57 a.m.

I've now finally finished my Polynomial RootSolver for the WP34s, and I think it's a quite good piece of software. :-)

Unlike almost every other such program (at least for calculators) it also calculates multiple roots with full precision, where the usual programs fail miserably and give completely wrong roots.

Although I've tested my program with MANY example polynomials I've found in the dozens of sources (about polynomial root solving) that I've read (and so I far I get excellent results!), but before releasing my program I would like to make a few more tests.

So does anyone know a special collection of 'hard to solve' polynomials, especially with multiple roots?
But please no polynomials with multiplicity=degree, i.e. P(x)=(x-x0)^n, because this is too easy and is automatically solved exactly by the used Laguerre method. ;-)

Franz

      
Re: Polynomial roots again ...
Message #2 Posted by Werner on 23 Dec 2011, 1:25 p.m.,
in response to message #1 by fhub

This PDF http://www-irma.u-strasbg.fr/~bugeaud/travaux/PolSurvRev1.pdf gives some classes of polynomials with integer coefficients, yet arbitrarily close roots.

Cheers, Werner

            
Re: Polynomial roots again ...
Message #3 Posted by fhub on 23 Dec 2011, 2:33 p.m.,
in response to message #2 by Werner

Quote:
This PDF http://www-irma.u-strasbg.fr/~bugeaud/travaux/PolSurvRev1.pdf gives some classes of polynomials with integer coefficients, yet arbitrarily close roots.
Thanks Werner, but I was rather expecting a ready-to-use list of polynomials known to be a hard task for rootsolvers (like Valentin's matrix examples).

It's already enough work to enter lots of (usually quite long) coefficients into the WP34s, but I'm absolutely not keen on studying a long abstract mathematical paper just to be able to build such polynomials by myself - for this task I could just use a CAS.

Franz

Edited: 23 Dec 2011, 2:35 p.m.

                  
Re: Polynomial roots again ...
Message #4 Posted by Werner on 27 Dec 2011, 2:11 a.m.,
in response to message #3 by fhub

There's no need to study anything.. if you cared to scan the document for 5 minutes, you would see that the polynomial

x^d - 2*(a*x - 1)^2 = 0
has two roots separated by a distance less than a^(-(d+2)/2) If you take a to be a convenient 10, and d=30, for instance, this means that
x^30 - 200*x^2 + 40*x - 2 = 0
has two roots that agree to 16 digits, yet are distinct:
0.09999999999999992929
0.10000000000000007071
By making the order of the first term larger, you can get the roots arbitrarily close, but still distinct. This means that your routine to eliminate multiplicity in roots will not apply here, but numerically they will be treated as multiple roots, resulting in slower convergence and loss of accuracy. The HP48, working with 15 internal digits, returns
.0999999635205
.100000036479
losing *nine* digits.

Cheers, Werner

                        
Re: Polynomial roots again ...
Message #5 Posted by fhub on 27 Dec 2011, 6:44 a.m.,
in response to message #4 by Werner

Quote:
x^30 - 200*x^2 + 40*x - 2 = 0
has two roots that agree to 16 digits, yet are distinct:
0.09999999999999992929
0.10000000000000007071
By making the order of the first term larger, you can get the roots arbitrarily close, but still distinct. This means that your routine to eliminate multiplicity in roots will not apply here, but numerically they will be treated as multiple roots, resulting in slower convergence and loss of accuracy. The HP48, working with 15 internal digits, returns
.0999999635205
.100000036479
losing *nine* digits.
Hi Werner,

I've now played a bit with your polynomial above, and of course 'distinct' roots which agree to 16 digits are just beyond the calculation accuracy of PRS (or should I say WP34s?).

In line 402 of PRS there's my 'multiple root threshold' (SDL 10) which I've chosen after many tests with and without multiple roots. If I change that line to "SDL 08" then I get in fact 2 distinct roots for your polynomial, but (like the HP48) with much less significant digits - a 'multiple' root of 0.1 (which PRS returns with its default SDL 10) is definitely much more accurate than the 2 roots given by the HP48 (or if I change SDL to 08).

So I would indeed say that my program returns a much better result for this example. But as I already mentioned - if you construct polynomials with distinct roots (or root clusters) only in the 16th (or higher) digit, then you shouldn't wonder that you don't get absolutely exact results.

After your example I'm thinking about making this 'threshold' value for multiple roots (SDL 10) an option that can be changed by the user (like the accuracy 1e-15 in R99)!? OTOH this would make the usage of the program more complicated - and which user would really know when (and to which value) change this parameter for multiplicity testing ... :-(

Franz

Edited: 27 Dec 2011, 11:17 a.m.

                  
Re: Polynomial roots again ...
Message #6 Posted by Werner on 27 Dec 2011, 4:29 a.m.,
in response to message #3 by fhub

A bit further on, you see the more general formula

 x^d - 2*(a*x-1)^(d-1) = 0
that has d-1 roots very close together. For a=d=10, that amounts to
 x^10 - 2*(10*x - 1)^9 = 0
 x^10 - 2000000000*x^9 + 1800000000*x^8 - 720000000*x^7
 + 168000000*x^6 - 25200000*x^5 + 2520000*x^4 - 168000*x^3
 + 7200*x^2 - 180*x + 2 = 0
Here, the WolframAlpha and the HP-48 results differ in the third digit already.

Cheers, Werner

      
Re: Polynomial roots again ...
Message #7 Posted by Marcus von Cube, Germany on 23 Dec 2011, 1:40 p.m.,
in response to message #1 by fhub

With the most recent emulator you can test the double precision functionality. In DBLON mode you loose some registers, of course. With the full setup, 50 double precision registers are available. The set of lettered registers is reduced to X, Y, Z, T, L and I.

There are most probably some rough edges. I'd like to here from you.

            
Re: Polynomial roots again ...
Message #8 Posted by fhub on 23 Dec 2011, 2:29 p.m.,
in response to message #7 by Marcus von Cube, Germany

Quote:
I'd like to here from you.
Well Marcus, I'm not yet very impressed with the latest improvements (I'm still using SVN 2030 before the new flash handling has been introduced).

After my program had reached the 510 step limit I tried loading and running it in the newer versions, just to see that this reduces the number of available registers, and I definitely need many/all of them. Since I want to solve also very high-degree polynomials (my program limit is now 39), I'm using all registers R80-R99 for my calculations, so reducing the registers for more program steps is not acceptable for my program. So I finally decided to switch back to the older version again and to remove an (almost never used) feature of my program again, and now my final version is 504 steps, which is enough for polynomials upto degree 39 (and runs even in version 2 of WP34s).

I also don't know yet what I should think about this new DBL precision mode - is this really necessary (i.e. does it really give any better results)???
At least for my polynomial solver I don't think so, because also without DBL precision I already get results accurate to all displayed 12 digits, and I can't imagine that this new mode would do it any better. And with even more register reducement it's definitely a no-go for me.

And finally I'm absolutely not happy with this new library management - it was so much easier with single wp34s-n.dat files for each program!

Franz

Edited: 23 Dec 2011, 2:37 p.m.

                  
Re: Polynomial roots again ...
Message #9 Posted by Marcus von Cube, Germany on 23 Dec 2011, 2:55 p.m.,
in response to message #8 by fhub

I have to disagree. The new library management is way better in getting the most out of the rare flash memory.

Many registers and long programs work well together if you let the program run off library or backup space and use RAM exclusively for data with a combination of global and local registers. Local registers are very handy for internal variables while global registers should be used for user data.

The reason to introduce double precision mode was an idea to implement more functions in user code which will save space, hopefully more than the new feature costs. Double precision comes in handy when an algorithm is susceptible to cancellation. The more digits the better. The internal precision is already at 39 digits and we won't loose too much of the potential when we start implementing stuff in user code with the help of double precision registers.

                        
Re: Polynomial roots again ...
Message #10 Posted by fhub on 23 Dec 2011, 3:31 p.m.,
in response to message #9 by Marcus von Cube, Germany

Yes, if the reason for double precision is turning built-in functions into external user code, you're right of course. But for usual programs I see not need for it.

And my problem with the new library is that I can't get the libmanager working here (with my TinyPerl), and so I still prefer separated dat-files created by the perfectly working assembler. But of course it's just a matter of taste, and maybe I'll like it more when I become familiar with it anytime in the future ...

Edited: 23 Dec 2011, 3:32 p.m.

      
Re: Polynomial roots again ...
Message #11 Posted by Valentin Albillo on 23 Dec 2011, 4:49 p.m.,
in response to message #1 by fhub

Quote:
I've now finally finished my Polynomial RootSolver for the WP34s, and I think it's a quite good piece of software. :-) [...] So does anyone know a special collection of 'hard to solve' polynomials, especially with multiple roots?

Last time I heard of your "Polynomial RootSolver for the WP34s" was about a month ago, in this post of yours (the highlighting is mine):

    "Message #24 Posted by fhub on 20 Nov 2011, 3:57 p.m.,
     in response to message #23 by Dieter

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

So I take it that now you're asking for people to share their knowledge and materials with you ("a special collection of 'hard to solve' polynomials") while publicly stating that you won't share a thing with anyone because of some bruised ego you've got or something ?

Or have you (silently) changed your mind since that post of yours ?

Not that a change of mind on your part would be unheard of, actually ...

Best regards from V.

            
Re: Polynomial roots again ...
Message #12 Posted by fhub on 23 Dec 2011, 5:48 p.m.,
in response to message #11 by Valentin Albillo

Would you really say that posting a few polynomials (or a link to them) is "share their knowledge and materials" with me?
Well, then I've also "shared my knowledge" more than a dozen of times in the past month when I reported bugs in the latest WP34s development. ;-)

And if you had read my posting a bit more mindfully then your reply won't have been necessary at all - I quote myself:
"..., but before releasing my program I would like to make a few more tests."

                  
Re: Polynomial roots again ...
Message #13 Posted by C.Ret on 28 Dec 2011, 12:42 p.m.,
in response to message #12 by fhub

Please forbid my poor english understanding. I am not native English reader. Can you explain me a bit; "releasing your program" does it mean that it will be soon possible to buy it?

All software people and companies I know speak or claim about “release” when they want to make me buy a new version of their products!

                        
Re: Polynomial roots again ...
Message #14 Posted by fhub on 28 Dec 2011, 12:49 p.m.,
in response to message #13 by C.Ret

Quote:
Please forbid my poor english understanding. I am not native English reader. Can you explain me a bit; "releasing your program" does it mean that it will be soon possible to buy it?

All software people and companies I know speak or claim about “release” when they want to make me buy a new version of their products!


LOL, of course not, I'm not a salesman or a company. ;-)
I'm also not a native English speaker, and maybe 'publish' would have been a better word than 'release'.

My program is already out (free of course) - see this thread:

http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/forum.cgi?read=207771#207771

Regards,
Franz

Edited: 28 Dec 2011, 12:49 p.m.


[ Return to Index | Top of Index ]

Go back to the main exhibit hall