Post Reply 
RPL second impressions (HP 28)
05-27-2022, 06:07 PM
Post: #64
RE: RPL second impressions (HP 28)
(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:
* ************************
* 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=-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=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):
        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):
        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)))
        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=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
* |B|>=.7
        GOSUB  EXab2     P
        ?ST=1  0         doing acos(z)?
        GOYES  docsc8
        GOSBVL =ACOS15   acos(P)        
        A=C    S         put sign(B)
        GOTO   docsc9
docsc8  GOSBVL =ASIN15   asin(P)
        ?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)
* -- done.
        GOSUBL rcscr
        GOSUBL stcd1     store imag part into (R0,R1)

(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.
Find all posts by this user
Quote this message in a reply
Post Reply 

Messages In This Thread
RPL second impressions (HP 28) - mdunn - 06-27-2018, 01:19 AM
RE: RPL second impressions (HP 28) - mdunn - 06-27-2018, 01:58 PM
RE: RPL second impressions (HP 28) - mdunn - 06-27-2018, 04:06 PM
RE: RPL second impressions (HP 28) - mdunn - 06-27-2018, 05:11 PM
RE: RPL second impressions (HP 28) - mdunn - 06-27-2018, 07:45 PM
RE: RPL second impressions (HP 28) - mdunn - 06-28-2018, 08:48 PM
RE: RPL second impressions (HP 28) - ttw - 07-04-2018, 10:52 PM
RE: RPL second impressions (HP 28) - Thomas Klemm - 05-27-2022 06:07 PM

User(s) browsing this thread: