How to implement trig-functions
*******************************
* ray presents *
* implementing trig-functions *
* may 2001 *
*******************************
*******************************
* ray presents *
* implementing trig-functions *
* may 2001 *
*******************************Intro
Hey fellas, time for a new tutorial, 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 within 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] which 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| < 1mirror 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| < 1quite 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
--------------------------------------------------------------------------------
