HP-71B Enhanced Math LEX

[J-F Garnier]

The HP-71B Math ROM is a great module, unfortunately it lacks some features that, surprisingly, are in its predecessor the HP-75C Math ROM, see here for a comparison.

A few months ago I started my "Math ROM source file" project to document the Math ROM code and I reached a point where I can start to implement enhancements. The result is a preliminary version with these additions:

- Support of complex arguments with the inverse trigonometric ASIN/ACOS/ATAN and hyperbolic ASINH/ACOSH/ATANH functions. For some unknown reasons these inverse functions didn't accept complex arguments in the HP-71B Math ROM. This is surprising because there is a lot of space available in the ROM (about 4KB) and the algorithms were well known by HP at the time since they are implemented in the HP-15C and in the HP-75C Math ROM too.
Using the HP-75C Math ROM source files as a guide to identify and document the algorithms, it has been surprisingly easy to implement the new functions in the HP-71B Math LEX after having the scalar real and complex routines properly documented.

- Addition of the cotangent secant and cosecant functions (COT SEC and CSC). These functions are present in the HP-75C (and in the previous series 80) and were initially included in the HP-71B mainframe. However they have been removed during the HP-71B development: a trace can be found in the HP-71B source files, AB&FCN module (Functions) line 1692, with this comment: "COT CSC & SEC removed from mainframe - 3/31/83 S.B.". Now they are back !

These new functionalities actually use three different enhancement techniques:
the complex ASIN/ACOS/ATAN use the poll mechanism to enhance mainframe keywords,
the complex ASINH/ACOSH/ATANH extend existing Math LEX keywords,
the COT, SEC and CSC functions are provided as new keywords in the Math LEX.
All these three techniques are now well documented and mastered.

I'm now thinking to add several matrix functions that are missing in the HP-71B Math ROM to reach at least the level of the HP-75C equivalent ROM. It will require understanding and documenting a lot of presently poorly annotated code.

You can find this preliminary test version 2a on the new "Math ROM page" on my site: http://www.jeffcalc.hp41.eu/emu71/mathrom.html
Feel free to test it with your favourite Math routines and of course your feedbacks are welcome.
I'm confident that I didn't break anything by making a new build of the Math LEX but still be careful and save your work first!

J-F

[rprosperi]

Wow, exciting announcement Jean-Francois, thank you for sharing this here!

I am tied up for a few weeks on another project, but look forward to checking this out in detail as soon as I can. Somewhere (!?) I have some very old notes on things missing/to be improved in the 71B Math ROM, I will look around for them. I believe the source of these notes was a PPC Computer Journal article and will look there as well.

Thanks!

[Valentin Albillo]

.
Hi, Jean-François:

(01-09-2020 07:30 PM)J-F Garnier Wrote:  A few months ago I started my "Math ROM source file" project to document the Math ROM code and I reached a point where I can start to implement enhancements. The result is a preliminary version with these additions: [...]

Wow, what a very nice surprise !!

As you pretty well know, I'm truly very fond of the HP-71B Math ROM and its awesome instruction set. I've been having a look at your new page dedicated to it and really you've made my day, it's an awesome project, a dream come true !

In the past I've tried to find source code for the Math ROM, pulling a few strings here and there, HP included, but to no avail as the source code seems to have been lost or else is buried under forty keys and impossible to get to.

My main interest was the source code for matrix operations, including documentation on parse routines, element handling via direct access, etc. I figured that if I could *see* such code I would be able to create my own specialized MAT keywords, as there's been many times where I had to use slow user-code loops (FOR..NEXT and other looping constructs) which could be replaced by a MAT keyword doing it all at assembly-language speeds.

I need to read your page more extensively and intend to see if I can contribute to it in any way. Nevertheless, I've enjoyed it immensely and of course want you to succeed in achieving the goals you state. Perhaps a list of proposed additional functionalities could be created to focus the implementation priorities.

Thank you very much for starting this wonderful project.

Very best regards.
V.

[charger73]

what a fantastic job... a great source of assembler coding...
thank you

(btw: who builds a sasm version for the hp75?)

[Raymond Del Tondo]

(01-10-2020 07:23 AM)charger73 Wrote:  (btw: who builds a sasm version for the hp75?)

AFAIK the 75C uses a CMOS version of the CPU used in the Series 80 desktops, so maybe the curator of Series80.org has info about suitable development systems?

[Sylvain Cote]

Thank you Jean-Francois for sharing your hard work and great findings. Impressive!
Sylvain

[Albert Chan]

Thank you for the enhanced Math Lex

This is a cubic solver, AX³ + BX² + CX + D = 0, using the new ASINH.
Code shift and scale to match identity: sinh(3z) = 3 sinh(z) + 4 sinh(z)^3

Code:

10 DESTROY ALL @ COMPLEX A,B,C,D,E,F,G,H,K
20 INPUT "A,B,C,D ? ";A,B,C,D
30 E=3*A @ H=-B/E
40 E=((E*H+2*B)*H+C)/A
50 F=(((A*H+B)*H+C)*H+D)/A
60 K=SQRT(E/.75)
70 E=ASINH(-4*F/K^3)/3 @ G=(0,2*PI/3) @ F=E+G @ G=E-G
80 E=SINH(E) @ F=SINH(F) @ G=SINH(G)
90 DISP E*K+H @ DISP F*K+H @ DISP G*K+H

>RUN
A,B,C,D ? 1,-6,11,-6
(2,0)
(1,0)
(3,-0)

>RUN
A,B,C,D ? 1,2,3,4
(-1.65062919144,0)
(-.174685404282,1.54686888723)
(-.174685404282,-1.54686888723)

[Albert Chan]

Bug report:

>SINH(ASINH((.5,0))
ERR:Invalid Expr

>SIN(ASIN((.5,0))
ERR:Invalid Expr

I am not sure if this is considered a bug ...

>DEGREES
>ASIN(.5)
30

>ASIN((.5,0))
(.523598775598,-0)
>SIN(RES)
(.5,-0)

[J-F Garnier]

(01-10-2020 05:01 PM)Albert Chan Wrote:  >SINH(ASINH((.5,0))
ERR:Invalid Expr
>SIN(ASIN((.5,0))
ERR:Invalid Expr

There is a missing closing ) in each of your expressions :-)

Quote:>DEGREES
>ASIN(.5)
30
>ASIN((.5,0))
(.523598775598,-0)
>SIN(RES)
(.5,-0)

Trig/hyperb/log/exp functions with complex arguments always use radians as unit.
This is consistent with other HP calculators/computers.
Only the ARG , POLAR and RECT functions are using the current angle mode.

Thanks for your feedback !

J-F

[Albert Chan]

(01-10-2020 05:40 PM)J-F Garnier Wrote:  Trig/hyperb/log/exp functions with complex arguments always use radians as unit.

Thanks for clearing this up. Sorry for the noise.

This equivalent cubic solver use identity: sin(3z) = 3 sin(z) - 4 sin(z)^3

Code:

10 DESTROY ALL @ COMPLEX A,B,C,D,E,F,G,H,K
20 INPUT "A,B,C,D ? ";A,B,C,D
30 E=3*A @ H=-B/E
40 E=((E*H+2*B)*H+C)/A
50 F=(((A*H+B)*H+C)*H+D)/A
60 K=SQRT(E/-.75)
70 E=ASIN(4*F/K^3)/3 @ G=2*PI/3 @ F=E+G @ G=E-G
80 E=SIN(E) @ F=SIN(F) @ G=SIN(G)
90 DISP E*K+H @ DISP F*K+H @ DISP G*K+H

[Erwin]

Hi,

don't know if this is helpful. I have a MATH ROM version from Rosenbaum - see my post with the link: MATH ROM 1x
to the old forum. I use this in my FRAM71. Great to see yours enhancements - very great job

regards
Erwin

[J-F Garnier]

(01-13-2020 06:10 PM)Erwin Wrote:  don't know if this is helpful. I have a MATH ROM version from Rosenbaum - see my post with the link: MATH ROM 1x to the old forum.

Yes, I know this version from Rodger Rosenbaum, since I was part of the 2007 discussion.
Actually, I tried to contact Rodger a few weeks ago, through his private email on the forum, to inform him of my work and ask if he would like to share the changes he made at the time, but didn't get any reply. He seems to visit the forum rarely now.
Does somebody have some news from Rodger, or know how to contact him?

J-F

[J-F Garnier]

An update of the "Math 2" LEX, version 2A is available.

Since no problem was found in the test version 2a, I promoted the version to the regular 2A level, including en passant a few additional enhancements:
- complex number support extended to the decimal log LGT function, to be in line with other calculators such as the 15C, 28C/S, 42S, and later models.
- complex number support also added to the alternate keywords ASN/ACS/ATN/LOG10.
- enhancement of the complex square operation z^2, now internally computed as z*z to provide the same functionality and accuracy for complex numbers than the x^2 key on other calculators (15C, 28C/S, 42S).

Mode details on the z^2 operation enhancement:
With the Math ROM 1A, an operation such as z^2 is internally computed as exp(2*ln(z)).
With complex numbers, this can lead to modest accuracy in some cases, such as:
>(1E-9,1)^2
(-1,2.00000769645E-9)
This is not a bug, other calculators give similar results (exactly the same for the 28C/S, 42S, etc) when using the y^x key.

Now, Math 2A handles the special case when the power is the real number 2 and gives:
>(1E-9,1)^2
(-1,2E-9)
like other calculators using the x^2 function.
Actually, it is still possible to use the general involution formula by forcing a complex power:
>(1E-9,1)^(2,0)
(-1,2.00000769645E-9)
Moreover, the special case is not based on a pattern detected during the expression parsing but at execution, so an expression like
X=2 @ (1E-9,1)^X is also handled.

The weakness of the z^2 operation in the HP71 Math ROM was mentioned in this article (last pages) of W. Kahan, now it is fixed!

J-F

[Albert Chan]

(01-26-2020 10:54 AM)J-F Garnier Wrote:  The weakness of the z^2 operation in the HP71 Math ROM was mentioned in this article (last pages) of W. Kahan, now it is fixed!

Thanks for the update !

Replacing Z^2 as Z*Z almost always give better result.
It take a while to find a counter-example.

>X=.254661579433 @ Y=.260740997689
>COMPLEX Z @ Z = (X,Y)
>Z^2
(-3.13334783654, .132801428589)
>Z*Z
(-3.13334783655, .132801428589)

Real part of Z = X*X - Y*Y = (X+Y) * (X-Y) gives better result.
> (X+Y) * (X-Y)         ! exact value = -3.13334783654 4934739232
-3.13334783654

[J-F Garnier]

(01-26-2020 03:08 PM)Albert Chan Wrote:  It take a while to find a counter-example.
>X=.254661579433 @ Y=.260740997689
>COMPLEX Z @ Z = (X,Y)
>Z^2
(-3.13334783654, .132801428589)
>Z*Z
(-3.13334783655, .132801428589)

Nice finding of a worst case for z*z, thanks! I believe that the z^2 result is better just by chance.
How did you find it? I'm guessing just brute force?
Anyway, the next 3 digits are 493, very close to the limit for up/down rounding, so we may say that both answers are equally right :-)

Quote:Real part of Z = X*X - Y*Y = (X+Y) * (X-Y) gives better result.
> (X+Y) * (X-Y)         ! exact value = -3.13334783654 4934739232
-3.13334783654

Interesting. The 71 is using the x*x-y*y formula, as do the following Saturn-based machines since they give the same answer for your worst case test.
I understand that the alternate formula (x+y)*(x-y) saves a multiplication, but does it guarantee that the result will be *always* better?

Thanks for your feedback!

J-F

[Albert Chan]

(01-27-2020 08:05 AM)J-F Garnier Wrote:  I understand that the alternate formula (x+y)*(x-y) saves a multiplication, but does it guarantee that the result will be *always* better?

Yes, I believe (x*x-y*y) is better calculated as (x+y)*(x-y)

|x| + |y| always accurate to whatever the system precision (unless it overflows)
|x| − |y| is exact if ratio of |x| and |y| = 1/2 to 2. (unless it underflows)
(see Kahans Miscalculating Area and Angles of a Needle-like Triangle, Section 4, Why Cancellation Cannot Hurt)

Below code searched bad REPT(Z*Z), for X = 0.1 to 1/3, 1 ≤ Y/X < 2

Code:

10 COMPLEX Z,Z1,Z2
20 X=.1+7/30*RND @ Y=X+RND*X @ Z=(X,Y)
30 Z1=Z*Z @ Z2=Z^(2,0) @ R=(X+Y)*(X-Y)
40 IF Z1=Z2 OR REPT(Z1)=R THEN 20
60 DISP "Z=";Z
70 DISP "Z*Z=";Z1
80 DISP "Z^2=";Z2
90 DISP "OK =";(R,2*X*Y)

>RANDOMIZE 1
>RUN
Z= (.235285090644,.236095767159)
Z*Z= (-3.82137391042E-4,.111099627953)
Z^2= (-3.82137391039E-4,.111099627953)
OK  = (-3.82137391041E-4,.111099627953)
>RUN
Z= (.180255457178,.190657390504)
Z*Z= (-3.85821071135E-3,6.87340701793E-2)
Z^2= (-3.85821071134E-3,6.87340701793E-2)
OK  = (-3.85821071134E-3,6.87340701793E-2)
>RUN
Z= (.1084548287,.127000096548)
Z*Z= (-4.36657465486E-3,.027547547432)
Z^2= (-4.36657465485E-3,.027547547432)
OK  = (-4.36657465485E-3,.027547547432)

Error (both Z*Z and Z^2) increases if |X| ≈ |Y|. After more runs, we have:

>SCI 12
>RUN
Z= (1.27043248466E-1 ,1.27045908048E-1)
Z*Z= (-6.75770947000E-7, 3.22806497255E-2)
Z^2= (-6.75770946606E-7, 3.22806497255E-2)
OK  = (-6.75770947060E-7, 3.22806497255E-2)

[J-F Garnier]

(01-27-2020 11:35 AM)Albert Chan Wrote:  

(01-27-2020 08:05 AM)J-F Garnier Wrote:  I understand that the alternate formula (x+y)*(x-y) saves a multiplication, but does it guarantee that the result will be *always* better?

Yes, I believe (x*x-y*y) is better calculated as (x+y)*(x-y)

|x| + |y| always accurate to whatever the system precision (unless it overflows)
|x| − |y| is exact if ratio of |x| and |y| = 1/2 to 2. (unless it underflows)
(see Kahans Miscalculating Area and Angles of a Needle-like Triangle, Section 4, Why Cancellation Cannot Hurt)

I tend to agree with you. (x+y)*(x-y) is better than (x*x-y*y) especially when x and y are close.
For instance:
>X=PI @ Y=PI-1E-11
>Z=(X,Y)
>Z*Z
(6.283E-11,19.7392088021)
>(X+Y)*(X-Y)
6.28318530717E-11
The last formula gives 8 more significant digits for the real part.
Note, however, than in this particular case, the involution formula gives an even worst result:
>Z^(2,0)
(6.27015817311E-11,19.7392088021)

I checked the x^2 (or SQ in RPL) function on HP Saturn-based machines up to the 49G+ and they all seem to use the x*x-y*y formula.
So I will refrain from implementing the (x+y)*(x-y) approach in the 71 Math 2 LEX, my goal is to align its performance with other machines of the time (75C, 28S, 42S), and the Z^2 case implemented as Z*Z already provides an interesting improvement.

J-F

[J-F Garnier]

(01-26-2020 10:54 AM)J-F Garnier Wrote:  An update of the "Math 2" LEX, version 2A is available.

There is a broken function in the January build with mixed real/complex multiplication:
>2*(3,4)
(6,0)
Sorry for that.

I fixed it and released an updated version, the February build identification is now:
>CAT MATHLEX
MATHLEX LEX 28527 02/08/20 11:34 3

J-F

[rprosperi]

Why not simply bump the version to 2b, and then 2B once confirmed? That's what revision letters are for and there are 24 more available.

[J-F Garnier]

(02-12-2020 01:55 PM)rprosperi Wrote:  Why not simply bump the version to 2b, and then 2B once confirmed? That's what revision letters are for and there are 24 more available.

A version 2b will come out soon, in a couple of weeks, with enhancements of the matrix function set.
I didn't want to spend a version letter for such an insignificant correction - nobody noticed it, right? :-)
And you know, I may need much more than 26 (or 52 with lowercase letters) versions if I start using one for each minor change, see how many builds were produced for Free42 or newRPL.

J-F