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

[Albert Chan]

It may help to explain PeterP's root squaring program.
This was a PM I sent to PeterP, and other members.

Albert Chan Wrote:Instead of min abs P roots, we solve for max abs Q root, Q(x) = P(1/x)

Graeffe's root squaring method (next row roots = previous roots^2)

We only show top 3 coefficients, because we only care for Q max abs root

q = [2,3,5,7,11,13,17,19,23,29, ...] // assumed infinite degree polynomial
→ [2^2, 11, 27]
→ [2^4, 95, 303]
→ [2^8, 671, 7775]
→ [2^16, 3.530559E6, 1.01291839E8]
→ [2^32, 8.11677068927E11, 4.099840909585279E15]
→ [2^64, 3.455854558672141E25, 1.816641529875401E31]
→ [2^128, −5.240706455638273E50, 2.748844184968018E62]
→ [2^256, 8.757340043013178E100, 7.559454552508339E124]
→ [2^512, 9.837400259693430E201, 5.714553957282261E249]
...

2nd column are not "pure squares" (note the negative sign), but 3rd column is.
Thus, roots to seek are complex conjugates.
Assuming roots are well separated now:

max abs Q root = surd(5.714553957282261E249/2^512, 512*2) = 1.2399046971021601
min abs P root = 1 / (max abs Q root) = 0.8065135992606103

---

[VA012b] is similar to [VA012a], a diffusion problem.
Root Squaring process also based from its neighbor cells.

b[n] = a[n]^2
b[n-1] = -a[n-1]^2 + 2*a[n]*a[n-2]
b[n-2] = +a[n-2]^2 - 2*a[n-1]*a[n-3] + 2*a[n]*a[n-4]
...

By the time big primes effect diffused to the top, its effects are miniscule.

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

[Albert Chan]

If we consider polynomial coefficients with geometric progression:

f(x) = 1 + (r*x) + (r*x)^2 + ... + (r*x)^n = ((r*x)^(n+1) - 1) / ((r*x) - 1)

From RHS numerator, all f roots abs = 1/r

P(x) = 2 + 3x + 5x^2 + 7x^3 + 11x^4 + ...

P roots, sorted in abs: (-0.6458±0.4832i), (0.4472±0.7248i), (-0.2853±0.8292i), ...
With corresponding abs: 0.8065, 0.8517, 0.8769, ...

Primes does not grow as fast as geometric progression. (sum of reciprocal primes diverges)

P min abs root is due to ratio, (5/2=2.5) > (11/5=2.2)
If the ratios were about the same, we expected f(x) roots pattern, with similar sized roots.

Example, R(x) = P(x)+0.5, 5/(2+0.5) = 2.

R roots, sorted in abs: (-0.6811±0.5122i), (0.4692±0.7114i), (-0.2905±0.8666i), ...
With corresponding abs: 0.8522, 0.8522, 0.9140, ...

If guess close enough, we can use Newton's method to zeroed in P min abs root.

[Valentin Albillo]

 
Hi, Albert Chan,

Thanks for your recent additional comments, I appreciate it. However, I have a thing or two to comment back, read on ... (all highlights are mine)

Albert Chan Wrote:

Valentin Albillo Wrote: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)

Slight error in logic.

Not at all, see below.

Albert Chan Wrote:If consecutive minimum abs both ~0.81, we cannot deduce trend apply to higher degree. (we can assume trend continues, but have to later test validity of assumption)

And I did test, it's just that I didn't want to make an already long post any longer by including unneeded data that most people won't be interested in, as they understand the scope of my articles and challenges and trust my results (which they can verify by themselves, if in doubt), but as you seem to like said data, here you are, the checks I did (which I saved but didn't post):

Degree  Min. Magnitude
-----------------------
  1     .66666 6666667
  2     .63245 5532034
  4     .65224 6975033
  8     .73550 1548954
 16     .76890 4440166
 32     .80252 6072477
 64     .80650 0035750
 96     .80651 362173
128     .80651 3599285

140     .80651 3599258
145     .80651 3599260
150     .80651 3599261
160     .80651 3599261
165     .80651 3599261
170     .80651 3599261
175     .80651 3599261
180     .80651 3599261
185     .80651 3599261
190     .80651 3599261 
195     .80651 3599261
200     .80651 3599261

and I was more than satisfied that the trend continued alright and converged to the correct solution, namely .806513599261.

Albert Chan Wrote:

Valentin Albillo Wrote: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.

Close, but not quite.

If this "minimum necessary degree" polynomial also gives abs ≤ 0.81, we are done.
However, if it gives bigger abs, we have to repeat again, with even higher degree polynomial.

But it didn't give "bigger abs", it gave  .806513599261, which is less than 0.81, so your comment doesn't apply at all and my statement is fully correct, not merely "close".

Now a word on the scope for my articles and challenges. In a past thread in which you took part (you posted 5 times no less) and so you surely read my posts there, I said:

"My articles are intended as just that, articles to be published in a physical fan-made magazine [...] or else on the WWW (MoHPC) for all kinds of fans, most of them not scholars, so my articles do not have the structure nor goals of a formal peer-reviewed paper."

In other words, my goal is first and foremost to provide entertainment to the forum members and HP calc fans in general, enticing them to think about the challenge and how to use their vintage HP calc to solve it, and if additionally they learn something new and interesting (and even useful) from my productions then so much the better. That is my goal.

Yours is obviously different, seemingly centered on lecturing, posting symbolic proofs and theoretical ramblings and lots of data obtained in Xcas sessions, lua, Mathematica, the works. Good for you and for the people who like (lots of) posts like that. Not my cup of tea here.

And please leave aside topics having little or nothing to do with my present Problem 2 (the references to integration, Borwein integrals and Kahan), save that for your own threads or where it's appropriate. Thanks.

V.