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

[C.Ret]

Ah! Ah!

It is with great pleasure and impatience that I discover this new stage of the challenge!

I rub shoulders with it, but I'm worried. The parity of the polynomial is a problem for me as well as the slowness of my HP-71B. Couldn't things be more complicated than they seem?

As I'm not very proud of my performance and the capabilities of my HP-71B (unless it's the other way around?) and in order not to waste too much time developing my solution, I decided to tackle a much simpler polynomial:

\(P_9(x)=2+3x+5x^2+7x^3+11x^4+13x^5+17x^6+19x^7+23x^8+29x^9\)

This polynomial has same resemblances with the polynomial P which interests us but has the enormous advantage of being infinitely shorter.

Please note that the derivate of \(P_9\) is easily determinate as:\(P_9'(x)=3+2\cdot 5x+3\cdot 7x^2+4\cdot 11x^3+5\cdot 13x^4+6\cdot 17x^5+7\cdot 19x^6+8\cdot 23x^7+9\cdot 29x^8\)

This allows me to code faster; Using Horner's scheme, Newton's root finding method and the prime number facility built in the JPC ROM cartridge, I manage to compose a code adapted to my HP-71B which in a few tens of seconds gives me the absolute value of the only real solution:


10 DESTROY ALL                 !! SRC#012b V.A. challenge (v.9)
 @ REAL D,E,K,P,X,Y
 @ X=-2/3 @ E=1.E-12             !  Initial guess and accuracy
20 REPEAT                     !! Computation of Y=P(X) and D=P'(X)
 @    P=29 @ Y=P @ D=0        !    P last prime     
30    FOR K=9 TO 1 STEP -1       !    Y/D computation loop
 @       D=D*X+K*P
 @       P=FPRIM(P-1,2)          !    FPRIM(P-1,2) return previous prime
 @       Y=Y*X+P 
 @    NEXT K                       
40    X=X-Y/D                 !  Newton's method X(n+1) = X(n) - P(x)/P'(x) 
 @ UNTIL ABS(Y/D)<E                             !    
 @ DISP ABS(X) @ END                            !! Display absolute value of root X


\( \left| x_9 \right| = .79462 \)

To increase speed and easiness, values of \( P_9 \) and its derivate \( P_9' \) are compute together in the same FOR TO NEXT loop.


My next step was to use this exact algorithm for the next polynomial:

\(P_{10}(x)=2+3x+5x^2+7x^3+11x^4+13x^5+17x^6+19x^7+23x^8+29x^9+31x^{10}\)


But I suddenly discover what already discover Werner and J.-F. Garnier; polynomials with only positive coefficients and even degree have no real root. So my previous version of the code didn't converge at all!

No stress, thanks to his Math module, this HP-71B is much more powerful that any basic calculator. It is very easy to modify the program to run the same algorithm but with complex values:


10 DESTROY ALL                 !! SRC#012b V.A. challenge (v.10)
 @ REAL E,K,P
 @ COMPLEX D,X,Y
 @ X=(.3,.7) @ E=1.E-12             !  Initial guess and accuracy
20 REPEAT                     !! Computation of Y=P(X) and D=P'(X)
 @    P=31 @ Y=P @ D=0        !    P last prime     
30    FOR K=10 TO 1 STEP -1      !    Y/D computation loop
 @       D=D*X+K*P
 @       P=FPRIM(P-1,2)          !    FPRIM(P-1,2) return previous prime
 @       Y=Y*X+P 
 @    NEXT K                       
40    X=X-Y/D                 !  Newton's method X(n+1) = X(n) - P(x)/P'(x) 
 @ UNTIL ABS(Y/D)<E                             !    
 @ DISP ABS(X) @ END                            !! Display absolute value of root X

For \( P_{10} \), I easily get the following absolute minimal value (in fact norm of the complex root).

\( \left| z_{10} \right| = .734576 \) since the closest complex root I found for \( P_{10} \) near \((0,0)\) is \( z_{10} = (0.297370,0.671694) \)


P.S.: Using this complex valued algorithm with \( P_{9} \) clearly indicate that I have not found the minimum absolute value among the roots since other complex roots exists that potentially have a smaller absolute norm.
For instance: \( \left| z_9 \right|= 0.71297 \) due to root \((-.56933,.42917)\) or conjugate.


Now, I need new batteries and a large bunch of time to run my code for the true challenge for the lengthy polynomial! 10'000 that's really huge!
I urgently need to find a way to determine the correct initial guess at first attempt.

References: Here is one of the documents that inspire me and help me develop this exact solution.
Newton’s method, and the fractal it creates that Newton knew nothing about. At 11:16 start the chapter about the 'Fun Facts', as fun as this Valentin Challenge. But you may watch there why I am stuck at the moment and need another good idea to efficiently start a 3-hour computation on the right initial guess!