A digression around VA's SRC #012b

[J-F Garnier]

In my first approach on the SRC #012b problem that Valentin posted recently, I attempted to find the real roots of the proposed 10000th degree polynomial:

(11-09-2022 04:05 PM)J-F Garnier Wrote:  There are (at least) two real roots !
There is a real root between -.996 and -.997 and another one between -.999 and -1, in accordance to my guess that real roots (if existing) would be close to -1.

Now that Valentin posted his original solution, I wanted to go back to the problem of finding the real roots.
I already located two roots, the question is now: is the HP-71B BASIC able to find the roots accurately and in a sensible time?
Contrary to the original problem and solution, there is no short cut, we have to actually compute all the 10000 terms of the polynomial, because the real roots are very close to -1.
The problem is finding an efficient and accurate way to evaluate the polynomial. Then, the solver of the Math ROM will do it.

Evaluating a 10000th degree polynomial.

There is, in principle, no difficulty with a piece of code such as:
DEF FNP1(X)
   N= 10000
   Y=2 @ P=1
   FOR I=1 TO N @ P=FPRIM(P+2) @ Y=Y+P*X^I @ NEXT I
   FNP=Y
END DEF

But this code provides moderate accuracy and runs quite slow.
I used several ways to optimize it:
- compute the terms by pairs (odd-even powers), since they are partially cancelling each other,
- use the DOT function (of the Math ROM) that computes a sum of products, exactly what we need here, with internal 15-digit accuracy,
- and finally precalculate the prime numbers, since several evaluation of the polynomial will be needed when searching for the real root(s).

At the end, I'm using about 30KB for the storage of the prime numbers, and 80KB for the two vectors needed for DOT.
This is the cost for the efficiency, but perfectly doable on a HP-71B or emulator with at least 128KB of RAM:

10  ! SRC12BJF
20  OPTION BASE 1
3   INTEGER D(10000) ! prime intervals, use about 30KB
40  DIM A(5001),B(5001) ! scratch arrays (about 80KB)
50  N=10000
60  !
70  ! setup: precalculate the prime intervals
80  DISP "SETUP..."; @ T=TIME
90  D(1)=1 @ P=3
100 FOR I=2 TO N @ Q=P @ P=FPRIM(P+2) @ D(I)=P-Q @ NEXT I
110 T=TIME-T @ DISP " DONE";T
120 !
130 ! find the real roots
140 T=TIME @ X1=FNROOT(-.996,-.997,FNP(FVAR)) @ T=TIME-T
150 DISP "X1=";X1;T @ DISP "----"
160 T=TIME @ X2=FNROOT(-.999,-1,FNP(FVAR)) @ T=TIME-T
170 DISP "X2=";X2;T @ DISP "----"
180 END
190 !
200 DEF FNP(X)
210   ! P,Q:primes, K,L: indexes, fill the arrays A() and B()
220   P=2 @ K=1 ! first prime
230   A(1)=2 @ B(1)=1 @ L=2 ! first coefficient
240   FOR I=1 TO N STEP 2
250     Q=P+D(K) @ P=Q+D(K+1) @ A(L)=Q+P*X @ K=K+2
260     B(L)=X^I
270     L=L+1
280   NEXT I
290   Y=DOT(A,B) ! evaluate the polynomial
300   ! DISP X;Y ! to trace the root search
310 FNP=Y @ END DEF


Here are the results obtained on Emu71/Win:

>RUN
SETUP... DONE 6.37 (seconds)
X1=-.996168277365 57.18 (seconds)
----
X2=-.999296380168 78.73 (seconds)
----
The setup uses about 6s, then the roots are found in about one minute each and are accurate to a few ULPs relative to the values I can get with Free42 using the polynomial evaluation code of Werner and the built-in solver.
Not bad for a modest 12-digit (15 internally) emulated machine.

A small digression inside the digression:
With his series of SRC, Valentin's intention is to show that "advanced vintage HP calcs [...] are NOW still perfectly capable of solving recently-proposed tricky problems".
On my side, I'm mainly interested to demonstrate that the Classic HP BASIC language is still an efficient and easy-to-use tool (and so easy to read or re-read), whatever the machine or emulator used to run it.
Computation speed and memory were quite limited on the HP-71B, but emulators running on modern computers greatly extent the usability of the HP-71B BASIC.

J-F