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 approximationIf 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 angleFlag 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 06LBL 05 # Handle alpha <= pi/6
XEQ 01 # Calculate S1 = sin(alpha)
STO* 07 # SinVal = SinVal * S1
GTO 08LBL 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 * S1LBL 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 00LBL 09
XEQ 01 # Calculate sin(beta)
STO* 07 # SinVal = SinVal * S2
GTO 00LBL 03 # Handle angle X <= pi/6
XEQ 01 # Call subroutine to calculate sin(x)
GTO 00LBL 04 # Handle X <= pi/2
3
/
XEQ 01 # Calculate Let S = sin(X/3)
XEQ 02 # Calculate S *(3 – 4 * S^2)
GTO 00LBL 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
*
RTNLBL 02 # Subroutine to calculate S * (3 – 4 *S^2)
STO 07
3
RCL 07
X^2
4
*
-
*
RTNLBL 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 angleFlag 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 6LBL 5 # Handle alpha <= pi/6
GSB 1 # Calculate S1 = sin(alpha)
STO* 7 # SinVal = SinVal * S1
GTO 8LBL 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 * S1LBL 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 0LBL 9
GSB 1 # Calculate sin(beta)
STO* 7 # SinVal = SinVal * S2
GTO 0LBL 3 # Handle angle X <= pi/6
GSB 1 # Call subroutine to calculate sin(x)
GTO 0LBL 4 # Handle X <= pi/2
3
/
GSB 1 # Calculate Let S = sin(X/3)
GSB 2 # Calculate S *(3 – 4 * S^2)
GTO 0LBL 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
*
RTNLBL 2 # Subroutine to calculate S * (3 – 4 *S^2)
STO 7
3
RCL 7
X^2
4
*
-
*
RTNLBL 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