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?

Taylor Series formulas:

Maclaurin Series formula:

arctan(x)=x - x^3/3 + x^5/5 - x^7/7 + x^9/9 - ... (-1 < x < 1)

\ 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 ;

\ 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)   ;