Series Approximations for Functions in
Forth
David L. Jaffe
|
|
In embedded environments it is sometimes necessary to
calculate a trig function. The options for doing this are:
- 1. use software floating point
- have to load most of package for just one
function
- have to convert to/from floating point
- 2. use hardware floating point
- have to load most of package for just one
function
- have to use a floating point processor
- have to convert to/from floating point
- 3. use a table
- may not be sufficient if precision is required
- fast execution
- 4. use series approximation
- how to implement?
- how many terms to calculate?
- what is the magnitude of error?
Trigonometry refresher
sinØ = Y/R
tanØ = Y/X
Ø = arctan(Y/X)
Series Formulas
Taylor Series formulas:
Maclaurin Series formula:
arctan(x)=x - x^3/3 + x^5/5 - x^7/7 + x^9/9 - ... (-1
< x < 1)
Check out arctangent series formula in
Excel
\ ARCTAN code - David L. Jaffe
\
\ arctan = (180/pi)[x - (x^3)/3 + (x^5)/5 - ... ] for -1 < x < 1 degrees
\ arctan = 90 - arctan(1/x) for -1 > x > 1 degrees
\
\ algorithm
\ pi = 355/113
\ x = tangent = n/d (n,d are integers)
\
\ scale = internal scale factor
\
\
\ 1st term = scale*(n/d)
\ 2nd term = (1st term)*(1/3)*(n/d)*(n/d)*(-1)
\ 3rd term = (2nd term)*(3/5)*(n/d)*(n/d)*(-1)
\ nth term = (n-1 term)*((2n-1)/(2n+1))*(n/d)*(n/d)*(-1)
\
: it ;
\
here " cls forget it load" constant reload
\
: load reload $>tib ;
\
0 value ntan \ numerator
0 value dtan \ denominator
0 value iter \ iteration counter
0 value arcsum \ accumulate result
28648 value scale \ 500 * (180)*(113/355)
1 value error \ run series until error <
-1 value itermax \ max iterations
\
: (arctan) ( n d - tenths-of-degree)
0 =: iter \ initialize iter#
2dup d0= if 2drop 0 \ arctan(0) = 0
else 2dup =
if 2drop 450 \ arctan(1) = 45 degrees
else =: dtan =: ntan \ initialize
scale ntan dtan */ dup =: arcsum \ first term
\ cr 0 . dup . cr \ display 1st term
begin incr> iter \ bump counter
iter 2* dup 1- swap 1+ */
ntan dtan */ ntan dtan */ negate
dup +!> arcsum \ add to total
\ iter . arcsum . dup . cr \ look at iter, sum, term
dup abs error < iter itermax = or \ error or iter limit?
until drop
arcsum 25 + 50 / \ round up and truncate
then
then ;
\
: tdisplay
iter 8 .r 7 .r cr ;
\
: arctantest cr cr cr
." Angle Iter Calc " cr
0 3 .r 0 1000 (arctan) tdisplay
5 3 .r 87 1000 (arctan) tdisplay
10 3 .r 176 1000 (arctan) tdisplay
20 3 .r 364 1000 (arctan) tdisplay
30 3 .r 577 1000 (arctan) tdisplay
40 3 .r 840 1000 (arctan) tdisplay
44 3 .r 966 1000 (arctan) tdisplay
45 3 .r 1000 1000 (arctan) tdisplay ;
\
\ n d abs(n)<abs(d) calc case
\ + + T 0 + arctan 0
\ + + F 90 - arctan 4
\ + - T 180 - arctan 2
\ + - F 90 + arctan 6
\ - + T 0 - arctan 1
\ - + F -90 + arctan 5
\ - - T -180 + arctan 3
\ - - F -90 - arctan 7
\
: arctan ( n d - tenths-of-degree)
2dup 2dup abs swap abs \ n d n d construct case#
< 1 and 2* 2* -rot \ n d f1 n d
0< 1 and 2* swap \ n d f1 f2 n
0< 1 and + + dup >r -rot r> \ f123 n d f123
3 > if swap \ transpose n d
then (arctan) abs swap \ arctan case#
case 0 of endof \ adjust for quadrants
1 of 0 swap - endof
2 of 1800 swap - endof
3 of -1800 + endof
4 of 900 swap - endof
5 of -900 + endof
6 of 900 + endof
7 of -900 swap - endof
." can't get here"
endcase ;
Results
Angle |
Iterations |
Calculated |
0 |
1 |
0 |
5 |
3 |
50 |
10 |
3 |
100 |
20 |
5 |
200 |
30 |
7 |
300 |
40 |
17 |
401 |
44 |
55 |
441 |
45 |
160 |
452 |
\ SINE code - David L. Jaffe
\
\ sin = [x - (x^3)/6 + (x^5)/120 - ... ]
\
\ algorithm
\ pi = 355/113
\ x = angle in radians = (degrees)*(355/113)*(1/180)
\
\ 1 degree = .0174532 radians
\
\ sin ranges from 0 to 1
\ output result * 10,000
\
\ degree scale factor = (10000/180)*(355/113) = 174.533
\
\ 1st term = x*(10000/180)*(355/113) = x*(500/9)*(355/113)
\ 2nd term = (1st term)*(1/6)*(x^2)*(-1)
\ 3rd term = (2nd term)*(1/20)*(x^2)*(-1)
\ nth term = (n-1 term)*(x/(2n+1))*(x/(2n))*(-1)
\
: it ;
\
: load " cls forget it load" drop 1- $>tib ;
\
0 value x
0 value iter \ iteration counter
0 value sinsum \ accumulate result
1 value error \ run series until error <
-1 value itermax \ max iterations
\
: (sin) ( 100ths of degree - 10000*sin ) \ -9000 to 9000
0 =: iter
5 9 */ 355 113 */ dup =: sinsum dup =: x
\ cr 0 . dup . cr \ display first term
begin incr> iter \ bump counter
x 10000 */
iter 2* dup rot + swap dup 1+ * / \ (term) / (2n) / (2n+1)
x 10000 */ negate
dup +!> sinsum \ add to total
\ iter . sinsum . dup . cr \ look at iter, sum, term
dup abs error < iter itermax = or \ error or iter limit
until drop sinsum ;
\
: sin ( 100ths of degress - 10000*sin ) \ -32767 to 32768
9000 /mod
case 0 of endof
1 of 9000 swap - endof
2 of negate endof
3 of negate endof
-1 of 9000 swap - negate endof
-2 of negate endof
-3 of 9000 swap - endof
endcase (sin) ;
|
|
|