Thank you Gerson Barbosa!!
Message #1 Posted by Namir on 24 Nov 2008, 9:47 a.m.
I was playing with the approximations of trigonometric functions and their inverses. The articles by Gerson Barbosa in this web site (especially the most recent one on implementing these functions on an HP-12C) proved to be most valuable. I had used many methods to calculate the tangent (some are based on calculating a good approximation for the sine function). The MinMax polynomial approximation that Gerson presents is definitely a winner!! I first implemented the algorithms in Excel VBA and then for the HP-41C and the HP-67. It all looks good!
Thanks Gerson for putting the time to share effective algorithms to calculate the sine and arc tangent functions.
Namir
Here first is the pseudo-code (all angles are in radians):
Given the angle X which is in the range of –pi to pi
Given the function BaseSinApprox(x) which is the MinMax polynomial approximation
If X = 0 then Return 0 as the sin(X)
If X < 0 then Chs = -1 else Chs =1
Let X = |X|
If X <= pi/6 then
Let SinVal = BaseSinApprox(X)
ElseIf X <= pi/2 then
Let S = BaseSinApprox (X/3)
Let SinVal = S * (3 - 4 * S^2)
ElseIf X <= pi then
Let alpha = X/2
Let beta = pi/2 - X/2
If alpha <= pi/6 then
Let S1 = BaseSinApprox(alpha)
Else
Let S = BaseSinApprox(alpha/3)
Let S1 = S * (3 - 4 * S^2)
End If
If beta <= pi/6 then
Let S2 = BaseSinApprox(beta)
Else
Let S = BaseSinApprox(beta /3)
Let S2 = S * (3 - 4 * S^2)
End if
Let SinVal = 2 * S1 * S2
End If
Return Chs* SinVal as sin(X)
Here is the HP-41C
Registered used:
R00 = angle X
R01 = angle alpha = X/2
R02 = angle beta = pi/2 - X/2
R03 = pi/6
R04 = pi/2
R05 = S
R06 = SinVal
R07 = Temp angle
Flag 0 = set when X < 0 and clear when X > 0
LBL "SINE"
LBL A
CF 00
"X?"
PROMPT
X=0?
GTO 00
X<0? # If X < 0 set flag 00
SF 00
ABS
STO 00 # Store |X|
PI
2
STO 07 # Initialize SinVal
/
STO 04 # Calculate and store pi / 2
3
/
STO 03 # Calculate and store pi / 6
X<>Y # Put angle X back in X register
X<=Y? # Is X <= pi/6?
GTO 03 # Calculate sin(X) using MinMax polynomial
RCL 04
X<>Y
X<=Y? # Is X <= pi/2?
GTO 04
2
/
STO 01 # Calculate and store alpha
RCL 04
X<>Y
-
STO 02 # Calculate and store beta
RCL 03
RCL 01
X<=Y? # is alpha <= pi/6
GTO 05
GTO 06
LBL 05 # Handle alpha <= pi/6
XEQ 01 # Calculate S1 = sin(alpha)
STO* 07 # SinVal = SinVal * S1
GTO 08
LBL 06 # Handle alpha > pi/6
3
/
XEQ 01 # Calculate S = sin(alpha/3)
XEQ 02 # Calculate S *(3-4 S^2)
STO* 07 # SinVal = SinVal * S1
LBL 08 # Handle beta
RCL 02
RCL 03
X<Y? # Is beta >= pi/6
GTO 09
X<>Y
3
/
XEQ 01 # Calculate S = sin(beta/3)
XEQ 02 # Calculate s * (3 – 4 * S^2)
STO* 07 # SinVal = SinVal * S2
GTO 00
LBL 09
XEQ 01 # Calculate sin(beta)
STO* 07 # SinVal = SinVal * S2
GTO 00
LBL 03 # Handle angle X <= pi/6
XEQ 01 # Call subroutine to calculate sin(x)
GTO 00
LBL 04 # Handle X <= pi/2
3
/
XEQ 01 # Calculate Let S = sin(X/3)
XEQ 02 # Calculate S *(3 – 4 * S^2)
GTO 00
LBL 01 # Subroutine to calculate sin(x)
STO 07
X^2
ENTER
ENTER
ENTER
ENTER
2.74201854577E-06
*
0.000198410347969
-
*
0.00833333320429
+
*
6
1/X
-
*
1
+
RCL 07
*
RTN
LBL 02 # Subroutine to calculate S * (3 – 4 *S^2)
STO 07
3
RCL 07
X^2
4
*
-
*
RTN
LBL 00 # Display the result
"SIN="
FS?C 00 # If flag 00 is set, change the sine
CHS
ARCL X
PROMPT
RTN
And Here is the listing for the HP-67:
R0 = X
R1 = alpha = X/2
R2 = beta = pi/2 - X/2
R3 = pi/6
R4 = pi/2
R5 = S
R6 = SinVal
R7 = Temp angle
Flag 0 = set when X < 0 and clear when X > 0
LBL A
CF 0
X=0?
GTO 0
X<0? # If X < 0 set flag 0
SF 0
ABS
STO 0 # Store |X|
PI
2
STO 7 # Initialize SinVal
/
STO 4 # Calculate and store pi / 2
3
/
STO 3 # Calculate and store pi / 6
X<>Y # Put angle X back in X register
X<=Y? # Is X <= pi/6?
GTO 3 # Calculate sin(X) using MinMax polynomial
RCL 4
X<>Y
X<=Y? # Is X <= pi/2?
GTO 4
2
/
STO 1 # Calculate and store alpha
RCL 4
X<>Y
-
STO 2 # Calculate and store beta
RCL 3
RCL 1
X<=Y? # is alpha <= pi/6
GTO 5
GTO 6
LBL 5 # Handle alpha <= pi/6
GSB 1 # Calculate S1 = sin(alpha)
STO* 7 # SinVal = SinVal * S1
GTO 8
LBL 6 # Handle alpha > pi/6
3
/
GSB 1 # Calculate S = sin(alpha/3)
GSB 2 # Calculate S *(3-4 S^2)
STO* 7 # SinVal = SinVal * S1
LBL 8 # Handle beta
RCL 3
RCL 2
X>Y? # Is beta >= pi/6
GTO 9
3
/
GSB 1 # Calculate S = sin(beta/3)
GSB 2 # Calculate s * (3 – 4 * S^2)
STO* 7 # SinVal = SinVal * S2
GTO 0
LBL 9
GSB 1 # Calculate sin(beta)
STO* 7 # SinVal = SinVal * S2
GTO 0
LBL 3 # Handle angle X <= pi/6
GSB 1 # Call subroutine to calculate sin(x)
GTO 0
LBL 4 # Handle X <= pi/2
3
/
GSB 1 # Calculate Let S = sin(X/3)
GSB 2 # Calculate S *(3 – 4 * S^2)
GTO 0
LBL 1 # Subroutine to calculate sin(x)
STO 7
X^2
ENTER
ENTER
ENTER
ENTER
2.742018546E-06
*
1.984103480E-4
-
*
8.333333204E-3
+
*
6
1/X
-
*
1
+
RCL 7
*
RTN
LBL 2 # Subroutine to calculate S * (3 – 4 *S^2)
STO 7
3
RCL 7
X^2
4
*
-
*
RTN
LBL 0 # Display the result
FS? 0 # If flag 0 is set, change the sine
CHS
CF 0
RTN
Edited: 26 Nov 2008, 7:51 p.m. after one or more responses were posted
| 
Go back to the main exhibit hall