RPL second impressions (HP 28)

[pier4r]

(07-04-2018 10:10 PM)Valentin Albillo Wrote:  That's the decades-long problem I've always resented: people will constantly produce this or that library, or this or that utilities ROM, or this or that utilities LEX files and so on and so forth, which are then published and/or made available for free or nearly so, intended to be used in people's own programs and such. So far so good.

Regrettably, most people either don't use them in actual, non-trivial programs or, if they do, they don't make them publicly available so in the end the community rarely gets full-fledged, useful programs but tons upon tons of routines, libraries, utilities, etc., apparently because people love to write libraries and utilities (it's like solving a puzzle or a challenge) but complex, useful programs not so much.

See for example the HP41 series' PPC ROM. Apart from the demonstration examples in its big manual I've very rarely seen "in print" any real application programs calling its functions as part of the program's code, nor have I ever called them from my 41C programs either.

Well writing libraries is always fun and they are used at least once (I hope), during the tests of those libraries. For the programs you are talking to, I am not sure I follow you.

While sharing a library or a function that can be applied to different problems is somehow useful, as it can be reused, sharing a program (or a solution) has a bit different scope. I like shared solution, even only to learn why this or that works how it works, but I can understand that people say "well, this solution may be interesting only to me, so I don't see the need to share it and make an effort in a post online". It is a pity, as someone else may find their effort instructive, but it happens often.

I, for one, share my solutions mostly on my open git repository that is not immediately reachable as it is one of the N git repos online. When I do some more interesting research I try to combine it in a sort of an article / post. Trying to cover the results (and on the way to get them) even just for myself of the future. As I will be able to use the results in a faster way if I organize them in an article.

Unfortunately more often than not I have little time so either I tinker a bit more or a write an article. Often I do the former instead of the more useful (as it pays off over time) latter part. Plus sometimes I want to read (for example the VA pdfs and what not) or do other tasks and so the time is gone.

[RMollov]

(07-04-2018 10:10 PM)Valentin Albillo Wrote:  The one and only RPL machine I own is precisely a new-in-box HP28S because I very much like the form factor and the keyboard. Alas, I've never used it and never will.

Sounds like a couple of oxymorons here.

[Thomas Klemm]

(07-01-2018 10:35 PM)Thomas Klemm Wrote:  That's great. Thus I assume these 109 lines of SysRPL code to calculate the complex inverse cosine are considered of "trivial size" as well.

I'm mostly interested in the HP-15C assembler code but I'm fine with the HP-71B as well.
Could you provide a similar analysis of this function so we could compare both results?

(02-28-2015 04:03 PM)Valentin Albillo Wrote:  The HP-71B Math ROM has been dumped for its use with several emulators/simulators, I have the dumped file itself.

In theory it could be disassembled, but without proper documentation it would be a gargantuan task to turn the resulting Saturn's assembler code into anything useful. Reverse-engineering the truly complex and highly-optimized algorithms of this particular 32K piece of code (50% of the size of the HP-71B ROMs themselves) would take lots of time to lots of skilled people and I don't think anyone would have that much time or interest for this ~30 year-old piece of code.

Well then, probably not.

[Thomas Klemm]

(01-09-2020 07:30 PM)J-F Garnier Wrote:  - 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.

This is the corresponding routine in the assembler code of release 2B7:

Code:

* ************************
* docasc (math2a) 
* computes asin(z) or acos(z)
* input:   z splitted in (A,B)=x,(R0,R1)=y
*          S0 clear if doing asin(z), set if acos(z)
* output: z in (A,B)=x, (R0,R1)=y
* ------------------------------
* algorithm (based on 75 sources, LG%COMp63 pdfp49):
*   Q = y^2 / ( sqr((|x|+1)^2+y^2) + (|x|+1) )
*   R = sqr((|x|-1)^2+y^2) +||x|-1|
*   S = y^2/R if R#0, 0 if R=0
*   M = Q+R if |x|>=1, Q+S if |x|<1    ; M is [2(A-1)] in 75's src
*   N = Q+S if |x|>=1, Q+R if |x|<1    ; N is [2(A-|x|)] in 75's src
*   imag part = sign(y,x)*log(1+M/2+sqr((M/2)*(M/2+2))
*             = sign(y,x)*log(1+(A-1)+sqr((A-1)*(A+1)) in 75's src
*      with sign(y,x)=sign(y) if y<>0, -sgn(x) if y=0
*      if doing acos(z), change sign of imag part
*   B = 2*x/(M+2)    ; =2x/[2A] in 75's src 
*   T = N/(M+2)      ; =[2(A-|x|)]/[2A] in 75's src
*   P = sqr(T*(2-T))     
* if doing asin(z):
*   real part = asin(B) if |B|<.7 , sign(B)*acos(P) if |B|>=.7
* if doing acos(z):
*   real part = acos(B) if |B|<.7 , asin(P) if |B|>=.7
*   if B<0 and |B|>=.7 then do real part = pi - real part
* ------------------------------
* W1, W2, W3 = scratch math stack levels
* M1, M2 = scratch registers (R0,R1), (R2,R3)
* ------------------------------
* J-F Garnier, dec.2019
* ************************
docasc  GOSUBL clrf      clear flags S1-S4
        ST=1   4         signal SB (inexact by default) for packc routine
        GOSUBL stscr             W1=x  save x for later use
        GOSUB  EXab1     AB=y    M1=x
* -- save sign(y,x) in R4[S]
        C=R0            
        C=-C-1 S         -sign of x 
        ?B=0   W         y=0?
        GOYES  docsc1      yes, done
        C=A    S           else get sign of y
docsc1  D=C    S
        C=R4
        C=D    S
        R4=C             R4[S]=sign(y,x)  
*        
* -- now computes Q
        GOSUBL square    AB=y^2        
        GOSUBL stab2             M2=y^2   
        GOSUB  EXab1     x        
        A=0    S         |x|
        GOSBVL =ADDONE   1+|x|
        GOSUBL stab1             M1=1+|x|
        GOSUBL square    (1+|x|)^2 
        GOSUBL rccd2
        GOSUBL ad15s     (1+|x|)^2+y^2
        GOSUBL sqr15     sqr((1+|x|)^2+y^2)
        GOSUBL rccd1
        GOSUBL ad15s     sqr((1+|x|)^2+y^2)+(1+|x|)
        GOSUBL rccd2
        GOSUBL xyex
        GOSUBL dv15s     Q=y^2/(sqr((1+|x|)^2+y^2)+(1+|x|))
*may check Inf/Inf xcptn here (due to y=Inf):
        ?XM=0
        GOYES  docscx
        GOSUBL rccd2
        GOSUBL xyex      Q=Inf
docscx  GOSUBL stscr     W1=Q  W2=x
*        
* -- now computes M, N
        GOSUBL rclw2     x
        A=0    S         |x|
        GOSBVL =SUBONE   do |x|-1 , not 1-|x|
        ST=1   5         S5=1 will indicate |x|>=1
        ?A=0   S         (|x|-1)>=0 ?
        GOYES  docsc2      yes, keep S5=1    
        ST=0   5           no, S5=0 iif (|x|-1)<0 i.e. |x|<1
docsc2  A=0    S         ||x|-1|=|1-|x||        
        GOSUBL stab1             M1=|1-|x||
        GOSUBL square    (|1-|x||)^2
        GOSUBL rccd2
        GOSUBL ad15s     y^2+(|1-|x||)^2
        GOSUBL sqr15     sqr(y^2+(|1-|x||)^2)
        GOSUBL rccd1
        GOSUBL ad15s     R=sqr(y^2+(|1-|x||)^2)+|1-|x||
        GOSUBL stab1             M1=R
        ?B=0   W         R=0?
        GOYES  docsc3      yes, skip y^2/R, use S=0 instead
        GOSUBL rccd2
        GOSUBL xyex
        GOSUBL dv15s     S=y^2/R  
*may check Inf/Inf xcptn here (due to y=Inf):
        ?XM=0
        GOYES  docsc3
        GOSUBL rccd2
        GOSUBL xyex      S=Inf
docsc3  GOSUBL rclw1     Q
        GOSUBL ad15s     Q+S
        GOSUBL stab2             M2=Q+S
        GOSUB  EXab1     R       M1=-
        GOSUBL rcscr     R Q     W1=x  
        GOSUBL ad15s     Q+R
        ?ST=1  5         |x|>=1?
        GOYES  docsc4
        GOSUB  EXab2     Q+S     M2=Q+R
docsc4  GOSUBL stab1             M1=M
* now AB=M1=M, M2=N
        GOSUB  get2.
        GOSUBL ad15s     M+2
        GOSUB  EXab2     N       M2=M+2
        GOSUBL stscr             W1=N W2=x
        GOSUB  EXab2     M+2 
        GOSUBL stscr             W1=M+2 W2=N W3=x   save for later use
*
* -- now computes imag part
        GOSUB  EXab1     M
        GOSUBL fdiv2     M/2
        GOSUBL stab1             M1=M/2
        GOSUB  get2.
        GOSUBL ad15s     M/2+2
        GOSUBL rccd1
        GOSUBL mp15s     (M/2)*(M/2+2)                
        GOSUBL sqr15     sqr((M/2)*(M/2+2))
        GOSUBL rccd1
        GOSUBL ad15s     M/2+sqr((M/2)*(M/2+2))
        GOSBVL =LN1+15   log(1+M/2+sqr((M/2)*(M/2+2)))
        C=R4
        A=C    S         put sign(y,x)
        ?ST=0  0
        GOYES  docsc5
        A=-A-1 S         negate for acos(z)
docsc5  GOSUBL stab2             M2=imag part (impt) ***
*
* -- now computes B and P
* here, W1=M+2 W2=N W3=x 
        GOSUBL rclw3     x
        GOSUBL fmult2    2*x
        GOSUBL rclw1  
        GOSUBL xyex      2*x  (M+2)
        GOSUBL dv15s     B=2*x/(M+2)
        GOSUBL stab1             M1=B **
        C=R4
        C=A    S
        R4=C             R4[S]=sign(B)
        GOSUBL rclw2     N
        GOSUBL rcscr     N (M+2)  
        GOSUBL dv15s     T=N/(M+2)
* save imag part from M2 to the math scr stack (before asin/acos)        
        GOSUB  EXab2     impt    M2=T
        GOSUBL stscr             W1=impt
        GOSUB  EXab2     T
        GOSUBL stab2             M2=T
        A=-A-1 S         -T
        GOSUB  get2.
        GOSUBL ad15s     2-T
        GOSUBL rccd2
        GOSUBL mp15s     T*(2-T)
        GOSUBL sqr15     P=sqr(T*(2-T))
        GOSUBL stab2             M2=P  **
*
* -- now computes the real part
        ST=1   =sRAD     set RAD mode for asin/acos
*                        Wrn: asin/acos use R0-R3, M1 and M2 not preserved        
        GOSUBL rccd1              
        GOSUBL xyex      B
        GOSUB  test.7    |B|<.7?
        GONC   docsc7          no, go next section
* |B|<.7 
        ?ST=1  0         doing acos(z)?
        GOYES  docsc6         
        GOSBVL =ASIN15   asin(B)
        GOTO   docsc9
docsc6  GOSBVL =ACOS15   acos(B)
        GOTO   docsc9
docsc7  
* |B|>=.7
        GOSUB  EXab2     P
        ?ST=1  0         doing acos(z)?
        GOYES  docsc8
        GOSBVL =ACOS15   acos(P)        
        C=R4
        A=C    S         put sign(B)
        GOTO   docsc9
docsc8  GOSBVL =ASIN15   asin(P)
        C=R4
        ?C=0   S         B>=0?
        GOYES  docsc9      yes, done
* doing acos(z), B<0 and |B|>=.7 then do rept=pi-asin(P)
        A=-A-1 S         -asin(P)
        GOSBVL =PI/2     
        D=D+D  W         get PI
        D=D-1  A         (better, see MAKEPI in SM&MTH)
        GOSUBL ad15s     PI-asin(P)
docsc9
* -- done.
        GOSUBL rcscr
        GOSUBL stcd1     store imag part into (R0,R1)
        RTN

(01-09-2020 10:36 PM)rprosperi Wrote:  Wow, exciting announcement Jean-Francois, thank you for sharing this here!

I couldn't agree more.

Thanks to the comments I can follow the source code even though I'm not familiar with the entry points.
Some of them can be guessed easily, like fdiv2, mp15s, ad15s, dv15s or sqr15.
The program is mostly a list of go-subroutine calls with a little bit of "real" assembler code sprinkled here and there.

To me that looks similar to the SysRPL code.
Entry points start with %% and instead of using assembler to deal with registers, it uses stack shuffling instructions.
Maybe not a surprise as both use the same algorithm.

I have to admit that I wasn't aware that the complex ACOS function was missing from the original HP-71B Math ROM.
This explains why I never got a response to my request.
I am very happy that I finally got it.

[rprosperi]

(05-27-2022 06:07 PM)Thomas Klemm Wrote:  [snip...]

To me that looks similar to the SysRPL code.
Entry points start with %% and instead of using assembler to deal with registers, it uses stack shuffling instructions.
Maybe not a surprise as both use the same algorithm.

It's not surprising since many many members that worked on the 71B OS and IDS went on to work on the 28C/S and then the 48 ROM. Most folks don't realize that a vast majority of the underlying math libraries/routines in all post 71B machines up through the 50g are based on the foundation laid by the 71B OS. Cyrille has even commented that much of the 39GII and Prime's "guts" re-uses much of the same stuff, though ported to C for portability.

[Didier Lachieze]

(05-27-2022 07:38 PM)rprosperi Wrote:  Most folks don't realize that a vast majority of the underlying math libraries/routines in all post 71B machines up through the 50g are based on the foundation laid by the 71B OS.

Bill Wickes provided a port of the HP 71 Math Pac PROOT function for the HP 48S/SX which is the same assembly language code with a RPL front end: HP 48 Polynomial Rootfinder

[J-F Garnier]

(05-27-2022 06:07 PM)Thomas Klemm Wrote:  The program is mostly a list of go-subroutine calls with a little bit of "real" assembler code sprinkled here and there.

To me that looks similar to the SysRPL code.
Entry points start with %% and instead of using assembler to deal with registers, it uses stack shuffling instructions.
Maybe not a surprise as both use the same algorithm.

Yes, it looks a bit similar but with important differences:
- RPL is using the stack for each operation.
- HP-71 assembly uses the CPU registers pairs (A,B), (C,D) for the 1 or 2 arguments of each operation, and two scratch registers formed with registers R0-R3.

The HP-71 also uses a 4-level stack for keeping intermediate values, but it's only for storage, operations are not directly done one the stack.
For more complex algorithms, as encountered in the matrix operations, the HP-71 also uses other temporary memory locations, that must be managed "by hand" by the programmer.
So programming in HP-71 assembly requires to juggle with these limited resources (with many opportunities for bugs), whereas System RPL, with its infinite stack, is much easier to program.

It was the main benefit (and objective I guess) of System RPL, but it has a cost: a significantly slower execution speed. It is particularly visible between the 32S and 42S, the 32S even with its slower 640 kHz CPU is faster than the 42S for many operations.
Although I recognize the benefits of System RPL vs assembly, I prefer to program in assembly for efficiency, in a way not so different from the many people who seem to still prefer programming in RPN/RPL rather than in a clear algebraic language :-)

Quote:I have to admit that I wasn't aware that the complex ACOS function was missing from the original HP-71B Math ROM.
This explains why I never got a response to my request.
I am very happy that I finally got it.

On my side, I missed your article, it would have been useful for me when writing the complex inverse trig functions for the HP-71 !

(05-28-2022 05:38 AM)Didier Lachieze Wrote:  Bill Wickes provided a port of the HP 71 Math Pac PROOT function for the HP 48S/SX which is the same assembly language code with a RPL front end: HP 48 Polynomial Rootfinder

Unfortunately, Bill just provided the binary code , not the source file. PROOT code is still not commented in my "HP-71 Math ROM source file project". It may not be too difficult because the method is properly documented in the HP75/71 Math ROM manuals, so I may do it someday when I will work again on the HP71 Math Pac.

J-F