home UCM22

-----

.tHE .sIRIUS .cYBERNETICS .cORPORATION


UCM 22
                        *******************************
                        *        ray presents         *
                        * implementing trig-functions *
                        *          may 2001           *
                        *******************************

Intro
-----

Hey fellas, time for a new tute, I guess...a short one this time. well this time
I'm gonna try  to cover a easy  way of implementing trigonometrical functions on
the M680x0. It might  be useful  for precalculating  lookuptables  and  stuff on
machines which lack a FPU whithin your own program.

But before, the  ordinary contact  lines (if you're  having problems, advices or
anything):

 email: reimund.dratwa@freenet.de
 hp:    http://www.ray-tscc.de.gs


Story
-----

Once upon a time...
I went into real  trouble  getting a good  offset-table  for a tunnel  which had
built at  runtime (the actual problem). This  problem  was based  on a  devilish
math-gnome called arctan (urghh, he looks like some toadstool with arms and legs
or something =) ) needed for that. Now some might come  up saying "Hey why din't
you end up just  including an atan table ?"...easy, the thing had to be as small
as possible.
Good, I know, there are  probably some  nasty tricks around to avoid usage of an
atan for a tube  lookup. But mine  should look  good, too and  I wanted it to be
mathematical  correct so I  chose the harder  way which didn't appear to be that
hard anymore later one.
I remembered a  way for complex math-functions approximations called the 'taylor
series'... using them you'll be able to approximate hyperbolic, trigonometrical,
inverse-trigonometrical functions and some more within a certain interval.
I'm just going to cover the common ones namely the trig and trig^-1 functions.

some math....

sin, cos and tan approximations:

 -> sin(x) = x - x^3/3! + x^5/5! - x^7/7! ... +/- x^n/n!   (n faculty)

 works fine for x = [-pi/2;+pi/2] wich is equal to an interval between
 -90° and + 90° (notice: x is arc not deg !)

 -> cos(x) = 1 - x^2/2! + x^4/4! - x^6/6! ... +/- x^n/n!

 same as above

 -> tan(x) = x + x^3/3 + 2x^5/15 + 17x^7/315 + 62x^9/2835 ...

             + [2^2n(2^2n-1)Bn x^2n-1]/(2n!)

 complicated, isn't it...but it works wiht |x| < pi/2

 -> cot(x) = 1/x - x/3 - x^3/45 - 2x^5/945 ...

             - [2^2n Bn x^2n-1]/(2n)!

 for x = [0;pi]


that on the trigs...and now for the inv-trigs:

 -> atan(x) = x - x^3/3 + x^5/5 - x^7/7 ...+/- x^n/n  |x| < 1

 mirror to pi/2 or take a look into the source if you wanna cover
 a full quadrant

 -> asin(x) = x + 1x^3/2*3 + 1*3x^5/2*4*5 + 1*3*5x^7/2*4*6*7... |x| < 1

 -> acos(x) = pi/2 - asin(x) |x| < 1

quite hard thing you might think...well, me too. but  I'm neither  gonna explain
how exactly they work nor trying to understand this right now :).
So you're not alone...


Implementation
--------------

Basically, these formulas are quite clear structured and thus easy to implement.
And the great advantage: they don't need any lookup-tables so things will remain
short, even reasonably  quick for  ST purposes (remember we do not have a FPU on
the ST, so it's better than nothing).
The more elements you're going to use the more accurate your results will be but
since we're dealing with a limited number of data-registers 5 to 7 elements have
to be enough if  we're not going  to use variables or stackspace, which would be
really getting slow, I think.
At the same hand  we have our  limited fixed-point  precision...so more elements
wouldn't make sense with that in mind.
I used 8.8 fixedpoint  precision for instance...with 6 elements this is gonna do
well for most of the functions.
For sin and cos you  might as well  use 2.14 fixedpoint, because results will be
in a range of -1;1...the more precision the better results.


Some code and closing words
---------------------------

That's all...keep  looking for  my new tutorials ;)...next  time maybe something
on raycasting a la wolf3d...bye.

For the last part I merged 2 examples on here to show  you how I implemented the
cos and the atan function, by now:

;--------------------------------- atan -----------------------------------

; *********************************************************
; * computes the arctan to a given tangens value          *
; *             - by ray//.tSCc. 2000 -                   *
; *********************************************************
; *                                                       *
; * input  D0.w  - tan in 8.8 fp                          *
; * output D0.w  - arc in 8.8 fp                          *
; *********************************************************

pi     EQU $0324  ; guess


arctan:        cmp.w  #256,D0   ; tan = 1
         beq.s  returnpi4
         bhi.s  adjtan     ; adjust tan
         bsr.s  precomp
         rts

adjtan:     move.l #$FFFF,D1
         divu   D0,D1
         move.w D1,D0      ; tan = 1/tan
         bsr.s  precomp

         move.w #pi/2,D1
         sub.w  D0,D1      ; arc = PI/2 - arctan(1/tan)
         move.w D1,D0      ; (mirrored)
         rts

returnpi4:    move.w #pi/4,D0    ; return PI/4
         rts


; *********************************************************
; * precomputes arctan for using a taylor-serie           *
; * at some places i had to cheat a bit to keep the error *
; * as small as possible :)                               *
; *********************************************************
; *                                                       *
; * input  D0.w  - tan in 8.8 fp                          *
; * output D0.w  - arc in 8.8 fp                          *
; *                                                       *
; * D1     tan^2                                          *
; * D2-D7 elements of taylor serie                        *
; *********************************************************

precomp:    movem.l D1-D7,-(SP)

         move.w D0,D1
         mulu   D1,D1
         lsr.w  #8,D1      ; D1 = tan^2

         move.w D1,D2
         mulu   D0,D2
         lsr.w  #8,D2
         divu   #3,D2      ; D2 = tan^3/3

         move.w D2,D3
         mulu   D1,D3
         lsr.w  #8,D3
         divu   #5,D3      ; D3 = tan^5/5

         move.w D3,D4
         mulu   D1,D4
         lsr.w  #8,D4
         divu   #7,D4      ; D4 = tan^7/7

         move.w D2,D5
         mulu   D1,D5
         lsr.w  #8,D5
         divu   #9,D5      ; D5 = tan^9/9

         move.w D3,D6
         mulu   D1,D6
         lsr.w  #8,D6
         divu   #11,D6     ; D6 = tan^11/11

         move.w D2,D7
         mulu   D1,D7
         lsr.w  #8,D7
         divu   #13,D7     ; D7 = tan^13/13

; D0 = tan - tan^3/3 + tan^5/5 - tan^7/7 + tan^9/9 - tan^11/11 + tan^13/14
; (taylor-series)

         sub.w  D2,D0
         add.w  D3,D0
         sub.w  D4,D0
         add.w  D5,D0
         sub.w  D6,D0
                        add.w  D7,D0

         ext.l  D0      ; to be sure :)
         movem.l (SP)+,D1-D7

         rts

;---------------------------------- cos -----------------------------------

; *************************************************
; * computes the cosine to x using a taylor-serie *
; *************************************************
; * input  - D0  8.8 fp                           *
; * output - D0  8.8 fp                           *
; *************************************************
; *      by ray//.tSCc. 2000                      *
; *************************************************

; cos a = 1 - a^2/2 + a^4/24 - a^6/720    ( a = [-PI/2;PI/2] )

cos:     move.w  D0,D1           ; a^2
                mulu    D1,D1
                lsr.l   #8,D1

                move.w  D1,D2           ; a^4
                mulu    D2,D2
                lsr.l   #8,D2

                move.w  D2,D3           ; a^6
                mulu    D1,D3
                lsr.l   #8,D3

                lsr.l   #1,D1           ; D1 = D0^2/2
                divu    #24,D2          ; D2 = D0^4/24
                divu    #720,D3         ; D3 = D0^6/720

                move.w  #256,D0
                sub.w   D1,D0
                add.w   D2,D0
                sub.w   D3,D0           ; D0 = 1 - D0^2/2 + D0^4/24 - D0^6/720

                add.w   D0,D0
                rts

eof
--------------------------------------------------------------------------------

UCM 22