[VA] SRC #012b - Then and Now: Root

[Valentin Albillo]

  
Hi, all,

Well, a full week has elapsed since I posted my OP, which has already passed the 1,300 views mark (and Problem 1 has exceeded 9,000 views already,) and I've got a number of solutions and/or comments, namely by Werner, Jean-François Garnier, C.Ret, PeterP, Fernando del Rey and Albert Chan. Thank you very much to all of you for your interest and valuable contributions.

Now I'll provide my original solution to this Problem 2 but first a couple' comments:

1) Some of you provided just code but no numerical results, others provided numerical results but no code, and most of you provided comments but didn't provide timings. As I said in my previous post, I'd consider proper etiquette and most useful for everyone if contributors would kindly post program code, a sample run, results and timings. Not a rule but it would help. And comments are always most welcome, of course.

2) Again, I'm mildly surprised that no one posted RPL solutions or RPL code of any kind. As for RPN code, PeterP made a most brave attempt to get a reasonably accurate result (as did Werner with Free42, though using capabilities not available on the vintage HP-42S,) but it might be that RPL people deem this challenge as too trivial for their powerful vintage RPL calcs. On the other hand, perhaps it might be that ... naw, never mind.

That said, these are my own approach and resulting original solution:

First of all, the main difficulty is the 10,000^th-degree. Were it a mere 100^th-degree polynomial, it would be quite trivial, but dealing with the full polynomial using PROOT, say, would require allocating 10,000 real elements for the coefficients (they don't fit as integers) plus another 10,000 complex elements for the roots, for a grand total of ~ 30,000 x 8 = 240 Kb, to which you must add the considerable memory that PROOT needs internally (21 x 10.000 + 261 ~ 210 Kb,) which adds up to ~ 450 Kb, surely exceeding maximum available RAM).

Then again, finding the 10,000 roots woud take ~ (10,000/100)²x the time required to deal with a 100^th-degree polynomial, which on a physical HP-71B is about 2,100 sec., so the big one would take 2,100 x 100² ~ 243 days. In other words, utterly unfeasible. So much for sheer brute force ...

Now, from theoretical considerations it's pretty obvious that none of the 10,000 roots can have an absolute value (aka modulus, aka magnitude) > 1, because then the highest term, 104,743 x^10,000, would dominate the sum, making it nonzero. Likewise, if the roots are too close to 0 then their powers will quickly tend to zero, thus failing to ever contribute enough negative value to compensate for the lowest coefficient, 2. Thus, all the roots must reside in an annulus of external radius 1 and internal radius to be determined by the minimum magnitude among all the roots, which we can roughly estimate as I'll explain in a moment.

Once we have a rough estimation for said minimum magnitude, we can compute another estimation, this time for the minimum degree of the truncated polynomial so that its highest term already contributes negligibly to the evaluation of the polynomial according to the precision required, say 12 digits.

My resulting program is thus this 5-liner: (REPEAT/UNTIL are from the JPC ROM, just for show)

1  DESTROY ALL @ @ OPTION BASE 0 @ M=0 @ W=0 @ N=0 @ FIX 2
2  REPEAT @ W=M @ N=N+10 @ GOSUB 4 @ DISP N;M;P @ UNTIL ABS(M-W)<.01 @ STD
3  N=IP((-12-LGT(P))/LGT(M)) @ DISP "Deg:";N;".." @ GOSUB 4 @ DISP "Min:";M;P @ END

4  DIM A(N) @ COMPLEX R(N) @ P=2 @ A(N)=P @ FOR I=1 TO N @ P=FPRIM(P+1) @ A(N-I)=P
5  NEXT I @ MAT R=PROOT(A) @ M=1 @ FOR I=0 TO N-1 @ M=MIN(M,ABS(R(I))) @ NEXT I @ RETURN

Line 1 does some initialization.

Line 2 finds the minimum absolute value for polynomials of degrees 10, 20, 30, ... until two consecutive minimum values are within 0.01, which essentially gives us the result correct to 2 digits (~0.81)

Line 3 now uses this 2-digit value to compute an estimation to the minimum necessary degree to get a fully accurate result, which it then obtains and displays, correct to 12 digits. It is important to use the minimum-degree truncated polynomial for speed reasons ("the faster, the better") because PROOT finds the roots of a 153^th-degree polynomial about 71% faster than for a 200^th-degree one, say.

Lines 4 and 5 are a subroutine which creates the N^th-degree polynomial, fills it up with the prime coefficients, computes all N complex roots and returns the minimum magnitude among them. All of it in just 2 lines of code.

This way, the user doesn't have to guess and provide any particular degree for the truncated polynomial at all, the program finds the necessary degree and then the sought-for minimum absolute value to maximum accuracy (12 digits). In other words, the program does all the work, fast.

Let's do a sample run:

>RUN

Deg   Min   Ncoef
-----------------
 10   0.73   31
 20   0.80   73
 30   0.81   127
 40   0.81   179
 Deg: 153 ..
 Min: .806513599261   887

correct to 12 digits in 1 hr 18 min on a physical HP-71B. We can check the highest term contribution like this:

>P;M;N;P*M^N  ->  887   .806513599261   153   4.56575832163E-12

which indeed contributes negligibly to the value of the polynomial so it and all other higher terms can be safely ignored.

As for using PROOT to find all the N roots (most of them complex) vs. using Newton's method to directly find the root that results in the desired minimum magnitude (as A.Chan does in post#18,) the reason is that in general that might not guarantee that you get the absolute minimum.

In A.Chan's post it works but in general it could be just luck, no guarantee at all without further considerations to try and find a working initial guess, and as you can see in the graph below most roots are within a thin annulus of radii 1 and 0.82, i.e. only 0.18 wide, and the conjugated pair we need are extremely near to it, so the roots are quite clustered and isolating the one we want might be a hit-or-miss affair, while using PROOT and then finding the minimum magnitude is 100% guaranteed to find the correct one, giving us perfect ease of mind with minimum effort on our part.



Well, that will be it for now. Thanks again to all of you who contributed and some of you even solved both Problem 1 and Problem 2, a perfect 2 for 2 score so far, but as they were the easiest problems of the lot I wonder if you'll manage to also solve incoming Problem 3, which is sure to raise the bar ... slightly ?  We'll see ...

V.