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

[Werner]

Why not, in FNP:

X2 = X*X

and in the loop:

B(L)=B(L-1)*X2

should be faster than 5000 exponentials?

Werner

[J-F Garnier]

(11-25-2022 02:21 PM)Werner Wrote:  Why not, in FNP:

X2 = X*X

and in the loop:

B(L)=B(L-1)*X2

should be faster than 5000 exponentials?

Yes, faster, but may introduce accumulated rounding errors, that are difficult to analyse in the final result. I chose the safest way.

J-F

[EdS2]

Thanks for digressing! Quite surprising to me that such a polynomial isn't just too unwieldy.

[Albert Chan]

(11-25-2022 02:00 PM)J-F Garnier Wrote:  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.

When I run your code in Emu71/DOS, I get hit with  "ERR L40:Insufficient Memory"
How do I add more memory to Emu71/DOS ?

Anyway, I removed DOT(A,B) code, and simply do horner's rule.
BTW, I wish we had HORNER, with the extra speed and precision.

I added two extra roots ±1, and solve below polynomial real roots instead.
Last term served only to straighten the curve, to speed up FNROOT a bit.

FNQ(x) = (2 + 3x + 5x^2 + 7x^3 + ... + p(n+1)*x^n) * (1-x) * (1+x)^3

10 DESTROY ALL @ OPTION BASE 1
20 N=10000 @ INTEGER D(N) ! PRIME INTERVALS, USE ABOUT 30kb
30 DISP "SETUP..."; @ T=TIME
40 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) ! = F(X)*(1-X)*(1+X)^3
70 Q=-P @ FOR I=N TO 1 STEP -1 @ Q=Q*X+D(I) @ NEXT I
80 Q=Q*X+2 @ FNQ=Q*(1+X)^3 @ END DEF
90 DISP "DONE",TIME-T
100 !
110 ! FIND THE REAL ROOTS
120 T=TIME @ X1=FNROOT(-.996,-.997,FNQ(FVAR))
130 DISP "X1=";X1,TIME-T
140 T=TIME @ X2=FNROOT(-.999,-.9999,FNQ(FVAR))
150 DISP "X2=";X2,TIME-T

>run
SETUP...DONE           102.86
X1=-.996168277368      177.95
X2=-.999296380168      208.74

[rprosperi]

(11-25-2022 11:14 PM)Albert Chan Wrote:  When I run your code in Emu71/DOS, I get hit with  "ERR L40:Insufficient Memory"
How do I add more memory to Emu71/DOS ?

Turn-off the emulated 71B (the device, not the emulator)

In the menu:

1. Click Edit... Port configuration
2. Click the drop-down to select a port (or just use Port1)
3. In module type drop-down, click RAM
4. In Size drop down, click 64K
5. Leave 'No. of chips' on Auto
6. Click Apply
7. Click OK

Turn-on the 71B and you should now have 64K more RAM than before. Depending on how you have EMU71 setup, you may have to confirm saving changes when exiting EMU71.

Use more calculators Albert! Emulators have their place, but for some things, there's nothing like a real machine.

[J-F Garnier]

(11-26-2022 03:18 AM)rprosperi Wrote:  

(11-25-2022 11:14 PM)Albert Chan Wrote:  When I run your code in Emu71/DOS, I get hit with  "ERR L40:Insufficient Memory"
How do I add more memory to Emu71/DOS ?

Turn-off the emulated 71B (the device, not the emulator)

In the menu:
...

Sorry, Bob, but Albert asked about Emu71/DOS

To add memory, edit the emu71.ini file and add the following lines in the [MODULE] section after the "0 RAM 32 ram00.bin" line:
0 RAM 32 ram01.bin
0 RAM 32 ram02.bin
0 RAM 32 ram03.bin
This will add 3x32k, which gives a total of 128K with the default 32k. Check with DESTROY ALL, MEM.
You can add more if needed.

J-F

[Albert Chan]

Hi, rprosperi, J-F Garnier

Thanks for the tips. With 128k, I am now able to run J-F code.

How do you figure out memory requirement ?
Example, how does INTEGER D(10000) translated to 30K ?
By trial and error, INETEGER type allowed -99999 to +99999

What about other declarations, say, SHORT, REAL, COMPLEX ?

[J-F Garnier]

(11-25-2022 11:14 PM)Albert Chan Wrote:  Anyway, I removed DOT(A,B) code, and simply do horner's rule.
BTW, I wish we had HORNER, with the extra speed and precision.

I added two extra roots ±1, and solve below polynomial real roots instead.
Last term served only to straighten the curve, to speed up FNROOT a bit.

FNQ(x) = (2 + 3x + 5x^2 + 7x^3 + ... + p(n+1)*x^n) * (1-x) * (1+x)^3

>run
SETUP...DONE           102.86
X1=-.996168277368      177.95
X2=-.999296380168      208.74

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.

Great trick to multiply the polynomial by 1-x, so only the differences of the prime numbers appear. It simplifies a lot.

J-F

[rprosperi]

(11-26-2022 10:20 AM)J-F Garnier Wrote:  Sorry, Bob, but Albert asked about Emu71/DOS

To add memory, edit the emu71.ini file and add the following lines in the [MODULE] section after the "0 RAM 32 ram00.bin" line:
0 RAM 32 ram01.bin
0 RAM 32 ram02.bin
0 RAM 32 ram03.bin
This will add 3x32k, which gives a total of 128K with the default 32k. Check with DESTROY ALL, MEM.
You can add more if needed.

J-F

Alas, the effort to try to assist is squelched by the inability to read...

Since Emu71/WIN is so much easier to run on modern Windows, I just assumed, rather than actually using my eyes... it can be such a serious limitation... <sigh>

Sorry for misleading you Albert, hope you didn't waste too much time struggling to figure out what the heck I was talking about. Menus??? Drop-down lists??? WTF???

[C.Ret]

(11-26-2022 03:18 AM)rprosperi Wrote:  Emulators have their place, but for some things, there's nothing like a real machine.

I completely agree, nothing like a real HP-71B. I turned mine on to time Albert's code, but since I'm a lazy user, I simplified his code a bit. I hope no one will blame me.

LIST
10 SETTIME 0 @ DESTROY ALL @ OPTION BASE 1 @ REAL I,N,P,Q,X @ N=10000 @ INTEGER D(N)
20 D(1)=1 @ P=3 @ FOR I=2 TO N @ Q=P @ P=FPRIM(P+2) @ D(I)=P-Q @ NEXT I
30 DEF FN FNF(X) @ Q=-P @ FOR I=N TO 1 STEP -1 @ Q=Q*X+D(I) @ NEXT I @ Q=Q*X+2 @ FNF=Q*(1+X)^3 @ END DEF
40 DISP P;TIME$ @ DISP FNROOT(-.997,-.996,FNF(FVAR));TIME$ @ DISP FNROOT(-.9999,-.999,FNF(FVAR));TIME$

RUN
104743 00:36:54
-.996168277368 01:48:48
-.999296380168 03:27:03


(This are times from one run only. Currently I am trying to run it again to check)

[J-F Garnier]

(11-26-2022 01:18 PM)rprosperi Wrote:  Since Emu71/WIN is so much easier to run on modern Windows, I just assumed, rather than actually using my eyes...

Actually, Emu71/Win is better suitable here for its higher speed.

(11-26-2022 11:24 AM)Albert Chan Wrote:  How do you figure out memory requirement ?
Example, how does INTEGER D(10000) translated to 30K ?
By trial and error, INETEGER type allowed -99999 to +99999

What about other declarations, say, SHORT, REAL, COMPLEX ?

From the manuals (!):
INTEGER range is +/-99999 and is coded on 3 bytes.
SHORT range is +/-9.9999E+/-499 and is coded on 4.5 bytes.
and REAL of course on 8 bytes.

COMPLEX and SHORT COMPLEX are just 2x the REAL/SHORT types.

J-F

[Albert Chan]

(11-25-2022 02:00 PM)J-F Garnier Wrote:  Here are the results obtained on Emu71/Win:

>RUN
SETUP... DONE 6.37 (seconds)
X1=-.996168277365 57.18 (seconds)
----
X2=-.999296380168 78.73 (seconds)
----

My ancient 22 years old Pentium 3 866MHz XP machine, running Emu71/DOS give similar performance.
(9.84, 60.23, 83.78)

In that sense, Emu71/Win is *very* slow.

I hope Emu71/DOS get update to use modern hardware, instead of running in DosBox, killing its performance.

[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

[C.Ret]

(11-26-2022 04:45 PM)Albert Chan Wrote:  I hope Emu71/DOS get update to use modern hardware, instead of running in DosBox, killing its performance.

I don't quite understand the point of further accelerating an emulator that is already over 110x the speed of real hardware.

Or, the goal is to demonstrate that the HP-71B is too old and too slow to be used for serious purposes now.

I strive to use only real material and I see that the slightest attempt takes several hours while others waste only a few minutes of their free time to obtain more results than I could ever obtain.

The HP-71B is like me, old and outdated.

---

EDIT 2022-11-27 09:43 (UTC+1 Paris):

Meanwhile, my second timing attempt is running, and I get the following times:

LIST
10 SETTIME 0 @ DESTROY ALL @ OPTION BASE 1 @ REAL N,I,P,Q,X @ N=10000 @ INTEGER D(N)
20 D(1)=1 @ P=3 @ FOR I=2 TO N @ Q=P @ P=FPRIM(P+2) @ D(I)=P-Q @ NEXT I
30 DEF FNF(X) @ Q=-P @ FOR I=N TO 1 STEP -1 @ Q=Q*X+D(I) @ NEXT I @ Q=Q*X+2 @ FNF=Q*(1+X)^3 @ END DEF
40 DISP P;TIME$ @ DISP FNROOT(-.997,-.996,FNF(FVAR));TIME$ @ DISP FNROOT(-.9999,-.999,FNF(FVAR));TIME$

DELAY 0,0
RUN
104743 00:36:54
-.996168277368 01:48:45
-.999296380168 03:12:26


These are hours, minutes and seconds. I hope this data will help you to adjust the calibrations of your emulators. And fortuitously that they will also help you to realize the energy and the patience which it is necessary to mobilize on the real material to try to follow your fentasques simulated fentasme.

---

EDIT 2022-11-27 10:26 (UTC+1 Paris):

(11-26-2022 08:06 PM)Albert Chan Wrote:  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.

Well this is what happens...

I barely have time to persist on my unique HP-71B. While it has not finished the first calculation, Albert Chan has already on his virtual DOS and WINDOWS development platforms had time to do many tests, research and progress...

..grrr..

I'll try his new algorithm next weekend, if I find the time.
I also stuck a post-it on the fridge's door to remind me to buy a new set of AAA batteries this Tuesday.
Albert Chan said it just took longer.

[J-F Garnier]

(11-26-2022 02:18 PM)C.Ret Wrote:  

(11-26-2022 03:18 AM)rprosperi Wrote:  Emulators have their place, but for some things, there's nothing like a real machine.

I completely agree, nothing like a real HP-71B. I turned mine on to time Albert's code, but since I'm a lazy user, I simplified his code a bit. I hope no one will blame me.
...
RUN
104743 00:36:54
-.996168277368 01:48:48
-.999296380168 03:27:03

Fantastic! I didn't expect so 'fast' results on the physical HP-71B.

(11-27-2022 02:36 AM)C.Ret Wrote:  

(11-26-2022 04:45 PM)Albert Chan Wrote:  I hope Emu71/DOS get update to use modern hardware, instead of running in DosBox, killing its performance.

I don't quite understand the point of further accelerating an emulator that is already over 110x the speed of real hardware.

Or, the goal is to demonstrate that the HP-71B is too old and too slow to be used for serious purposes now.

I strive to use only real material and I see that the slightest attempt takes several hours while others waste only a few minutes of their free time to obtain more results than I could ever obtain.

The HP-71B is like me, old and outdated.

As many of us, don't worry.

We should not oppose real 'old' machines to emulators.
As I explained a few times, for me the essence of a machine is in its firmware, so in this perspective an emulator is not different from the real machine.
See for instance the HP-15C LE: is it a 'true' HP-15C or an emulator? Actually both IMHO.

Matter of fact, my 'old' emulators don't attempt to reproduce all the subtitle hardware details, or even have the same look, but they are compatible enough to run the same code.

(11-27-2022 02:36 AM)C.Ret Wrote:  I also stuck a post-it on the fridge's door to remind me to buy a new set of AAA batteries this Tuesday.

Do you know that a wall adapter exists ? :-)
I can send you a spare one to support your efforts!

J-F

[C.Ret]

(11-27-2022 09:30 AM)J-F Garnier Wrote:  Do you know that a wall adapter exists ? :-)
I can send you a spare one to support your efforts!

Yes, it was while browsing the user manual that I discovered that there was an AC adapter.

This is a very nice proposal.

But I think it will not be useful, I have my HP-71B for two years now and I had never needed to change the batteries which were those installed by the previous owner.

No, my next investment for my HP-71B will be to take advantage of its HP-IL module to connect it to today's machines.

[Albert Chan]

Hi, C.Ret

Below give some idea of time it take, for a physical HP71B.

I am basing numbers from my "fastest" Pentium 3 866Mhz, 512MB XP machine, running Emu71/DOS
If Emu71/DOS can run on modern machine without DOS emulator, I expected another 10X speedup.

FNQ(X) version: 9.84 + 23.17 + 26.72 = 59.13

From your 3:12:26 and 3:27:03 timings, Emulator speed at 195X to 210X. Let's say 200X

Y=X+1 version: 9.84 + 50.51 + 58.00 = 118.35 --> HP71B take 6½ hours.

We could speed up a bit, assuming Q = IP(Q) + FP(Q) = q + r

70 FOR I=N TO 0 STEP -1 @ R=(Q+R)*Y-R @ Q=D(I)-Q+IP(R) @ R=FP(R) @ NEXT I

IP-FP version: 9.85 + 45.31 + 51.87 = 107.03 --> HP71B take 6 hours.

All 3 versions produce the same real roots.

---

My Toshiba laptop, Intel i3 M370 2.40 GHz 4 GB, running original FNQ(x) code.
DosBox SVN-lfn can pushed up to 105% (weird, but it does)

DosBox, max 50%: 102.86 + 177.95 + 208.74 = 489.55 --> 25X (HP71B speed)
DosBox, max 105%: 65.05 + 111.32 + 130.74 = 307.11 --> 40X (HP71B speed)

For other programs, my laptop run at 8X my old Pentium Box.
But laptop at max speed, Emu71/DOS + DosBox still run 5 times slower.

--> Cost of DOS emulation layer ≈ 8*5 = 40X

This matched J-F Garnier numbers.

(06-04-2020 08:18 AM)J-F Garnier Wrote:  As a comparison, Emu71/DOS runs at about x300 speed natively on my W2K 32-bit system, and at x8 speed in DOSbox on my W10 64-bit system, so a ratio of 40 due to the DOSbox emulation layer.

[rprosperi]

(11-27-2022 02:36 AM)C.Ret Wrote:  I strive to use only real material and I see that the slightest attempt takes several hours while others waste only a few minutes of their free time to obtain more results than I could ever obtain.

The HP-71B is like me, old and outdated.

Well, I can think of another similarity - both are appreciated. I enjoy your posts and learn much from them, please continue your outdated contributions...

[Albert Chan]

I guess the reason simple horner's rule work well is we really don't need super accurate y=f(x)
FNROOT is interpolating the secant line for root, from points (x1,y1), (x2,y2)

x = x2 - (x2-x1) * [(y2-0)/(y2-y1)]

As |x2-x1| get smaller, only a rough final ratio suffice.

Another reason is FNROOT, unlike secant's method, keep an "anchor" from the other side.
Convergence is slower, but avoid the problem of huge errors when f(x) closing to root.

Because of FNROOT, both FNQ(x) and FNQ(y=x+1) versions converge exactly the same path.
FNQ(y=x+1) improved accuracy is wasted.

   X                    FNQ(X)              FNQ(Y=X+1)           EXACT FNQ(X)
-.997             -3.21319628854E-7     -3.21319628895E-7     -3.21319628893E-7
-.996              1.12009058650E-7      1.12009058650E-7      1.12009058650E-7
-.9965            -1.72129040193E-7     -1.72129040185E-7     -1.72129040185E-7
-.996197103203    -1.74709075393E-8     -1.74709075547E-8     -1.74709075547E-8
-.9961705078      -1.36946078247E-9     -1.36946080158E-9     -1.36946080144E-9
-.996168448294    -1.05051545133E-10    -1.05051550303E-10    -1.05051550266E-10
-.996084224147     5.37864589010E-8      5.37864589106E-8      5.37864589104E-8
-.996168284115    -4.14702160931E-12    -4.14700483334E-12    -4.14700478100E-12
-.996126254131     2.63582466399E-8      2.63582466500E-8      2.63582466496E-8
-.996168277503    -8.29143191590E-14    -8.29276460508E-14    -8.29276562553E-14
-.996147265817     1.30466902823E-8      1.30466902845E-8      1.30466902845E-8
-.996168277369    -5.55826356659E-16    -5.64092304694E-16    -5.64057294087E-16
-.996157771593     6.49034644324E-9      6.49034650739E-9      6.49034650714E-9
-.996168277368     5.28822647438E-17     5.03647313957E-17     5.05964697450E-17