Thank you Gerson Barbosa!!



#2

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


#3

Very nice.

The pseudo code confused me at first, until I realized that most of your tests were meant to be "<=" not just "=".

Then the range reduction made sense. Looks like a neat job of RPN coding also.


#4

Sorry for the confusion.I have edited the original post to correct the errors (I was on my way out on a trip to another state for Thanksgiving when I posted my comments).

Namir

Edited: 24 Nov 2008, 11:19 p.m.

#5

You're welcome, Namir!!

But don't forget to take a look at the references I have listed in the 12C Platinum article. Without them I coudn't have made it!

Perhaps you'd like to take a look at the TurboBCD version:

http://www.geocities.com/gwbarbosa/prgms_4.html

Regards,

Gerson.


#6

Thanks for the Turbo pascal code!

Namir


Forum Jump: