(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.
Search
Member List
Calendar
Help
/
RPL second impressions (HP 28)
