A digression around VA's SRC #012b

[Albert Chan]

(11-26-2022 11:29 AM)J-F Garnier Wrote:  I don't like Horner method when dealing with thousands of terms, due to the possible accumulated rounding errors.
But it seems it works fine to find the roots here.

To improve accuracy, we split Q=q+r, and let y=1+X > 0

Q = Q*X + D[I]
(q+r) = (q+r)*(y-1) + D[I] = (D[i]-q) + ((q+r)*y-r)

We then force q an integer, fractional |r| ≤ 0.5

10 DESTROY ALL @ OPTION BASE 0
20 N=10000 @ INTEGER D(N) ! PRIME INTERVALS, USE ABOUT 30kb
30 DISP "SETUP..."; @ T=TIME
40 D(0)=2 @ D(1)=1 @ P=3
50 FOR I=2 TO N @ Q=P @ P=FPRIM(P+2) @ D(I)=P-Q @ NEXT I
60 DEF FNQ(X) @ Y=1+X @ Q=-P @ R=0
70 FOR I=N TO 0 STEP -1 @ R=(Q+R)*Y-R @ Z=IROUND(R) @ R=R-Z @ Q=D(I)-Q+Z @ NEXT I
80 FNQ=(Q+R)*Y^3 @ END DEF
90 DISP "DONE",TIME-T
100 !
110 ! FIND THE REAL ROOTS
120 T=TIME @ DISP "X1=";FNROOT(-.996,-.997,FNQ(FVAR)),TIME-T
130 T=TIME @ DISP "X2=";FNROOT(-.999,-.9999,FNQ(FVAR)),TIME-T

With 0.0007 < y < 0.004, this make horner's rule close to 15-digits precision.
We still ends up to the same roots, it just take longer.

X1 = -.996168277368
X2 = -.999296380168