How to calculate sinewave values in 680x0
It was somewhere when I was updating some 3d-routines that I realized: Blah! These old sinewave-fetchers just won't do anymore and I decided to standardize the whole lot into a bunch of macros. The goal was to make a size-optimised sinewave generator and a speed/size optimised sin/cos-fetcher.
For those of you that don't understand what the hell I'm on about: sinewave calculation is a thing very common in most demo-effects and games. Sinewaves can be used to create sound waves in samples, rotate 3d-point around axi, make pretty patterns for your sprites to move in, etc, etc.
In this article we will first look at the math and the basic algorithmic shit, before converting the algorithms to 68020 code.
Now there are two approaches to using sine values in your program:
- You calculate the sine-value every time.
- At the begin of your program you calculate an array full of consecutive sine-values and later on you can fetch the sine-values from this array in one or two instructions.
Because a the 68000-series up to the 68030 aren't equipped with hardware built- in sinewave-calculation you have to do it with a combination of 680x0- instructions. This is of course a very expensive operation in terms of processor time. So without an FPU or 040/060 CPU it would be wise to consider the second option! When no speed is needed the first one is ok, but when doing realtime stuff > please forget about the first option!!
How to do this? Well.. We check how sine-values are calculated by pocket- calculators :)
~~~~~~~~~~~~~~~~~~~~~~~~~~ :Math alert 1: ~~~~~~~~~~~~~~~~~~~~~~~~~
Any sinevalue can be generated by this formula:
x^3 x^5 x^7 x^9
sin(x) = x - --- + --- - --- + --- - ... + ...
3! 5! 7! 9!
Where "3!" denotes 3*2*1, "5!" denotes 5*4*3*2*1 and so on..
The "^" means "to the power of" in case someone didn't know already.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~In GFA-basic this formula can be easily explained. Just type in this in your interpreter and you will see even a function with only 4 terms can approximate the sine-value very precisely.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ :listing 1: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
x=0.2
PRINT "x := ";x
PRINT "OWN: ";x-x^3/(3*2)+x^5/(5*4*3*2)-x^7/(7*6*5*4*3*2)+x^9/(9*8*7*6*5*4*3*2)
PRINT "GFA: ";SIN(x)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~Note: as your "x" gets bigger, the accuracy of the sine-value decreases. You can see this in the GFA-program when using big values for x and comparing the "OWN" result to the "GFA" result.
So when you create a sinewave in an array, you should make sure you only calculate the first quarter of the wave with the formula. You can easily derive the rest of the sinewave from the previously calculated sinewave. You simply copy (and possibly negate) all the values into the next quadrants like so:
Quadrant I is y-mirrored into quadrant II and x-mirrored into quadrant III.
Quadrant IV is done by mirroring quadrant I in both axis. The picture says it all (I hope..).
ÿÿi sine.pi1 m 0 10 0 101ÿThis algorithm is used by one particular 68020-routine from the famous Seboz/Spirits. I've seen a lot of ATARI coders use it. And why not?
The sinewave-generation is quite efficient (only little instructions). And the fetching of a sine/cosine-value is a piece of cake.
I think he also does the mirroring stuff, but I'm not sure.
Furthermore I warn you: The code is very messy (sorry Seboz).
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ :listing 2: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
sin_gen:
lea sin,A0
moveq #0,D6
lea 512*2+4(A0),A1
movea.l A0,A2
moveq #512/4-1,D7
sin_gen1:
lea sin_fakdiv+16,A3
move.l D6,D5
moveq #7,D2
sin_gen2:
move.w D2,D1
move.l D6,D3
sin_gen3:
mulu.l D6,D4:D3
rol.l #8,D3
rol.l #8,D4
move.b D3,D4
move.l D4,D3
dbra D1,sin_gen3
divs.l -(A3),D3
add.l D3,D5
subq.w #2,D2
bcc.s sin_gen2
sin_gen4:
asr.l #1,D5
move.l D5,(A0)+
move.l D5,-(A1)
add.l #$06487ED5/512,D6 ;2*pi*(2^24)
dbra D7,sin_gen1
move.l #$007FFFFF,-(A1)
lea 512(A1),A1
move.w #512/2+512/4-1,D1
sin_gen6:
move.l (A2)+,(A1)
neg.l (A1)+
dbra D1,sin_gen6
rts
sin_fakdiv:
DC.L -6,120,-5040,362880 ;-3!, 5!, -7!, 9!
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~As you can see it doesn't require too much instructions and you've built up a complete sinewave-array with 512 elements.
Now all you have to do fetch a sinevalue is ...
* d0.w: angle (index to element in sinetable)
move.l (sin,d0.w*4),d1 * Fetch value in d1.lNote: if your angle (index) becomes bigger than 512 or less than 0, a modulo operation should be performed upon it. Otherwise the index would point to values beyond the sinewave!!
In assembler a modulo operation goes like this:
divu.w #512,d0
swap d0Or even better an optimised version. Works only for 2-1,4-1,8-1,16-1,32-1..
andi.w #512-1,d0All quite simple, but how do we fetch a cosine-value. This is very important for calculating sprite trajectories or doing point-rotations.
The theory is that you can fetch every cosine value from the same sine-wave! A cosinewave is basicly the same as a sinewave, but only shifted in phase 1/2pi (1/4 of complete tablelength).
So what do we do with these facts? Simple: we shift the phase of our index.
Like so:
addi.w #512/4,d0 * Shift phase by 1/2pi.
andi.w #512-1,d0 * Perform modulo.
move.l (sin,d0.w*4),d1 * Fetch value in d1.lBut hhhmmm.. This is two extra instructions..
EXACTLY. I faced the same problem some time ago and that's when I decided to make my own sinewave routines.
The ideas were to:
- Get rid of the extra instruction when fetching a cosine-value.
- Make a sine-generator that was even more efficient as Seboz's one.
- Put everything in MACRO so that the instructions could be easily used in my code.
- And possibly make it possible to fetch the sine and cosine in one go!
The first thing I did was getting rid of the longwords and instead I used normal words, because they are accurate enough anyway. I looked at some other (older) sine-generators and found out they used a special mechanism to enable faster fetching of cosine-values! They simply put the cosine-value next to the sine value in the table.
Instead of: Do this:
+----------------+ +----------------+----------------+
| sin(2pi/512*1) | | sin(2pi/512*1) | cos(2pi/512*1) |
| sin(2pi/512*2) | | sin(2pi/512*2) | cos(2pi/512*2) |
| sin(2pi/512*3) | | sin(2pi/512*3) | cos(2pi/512*3) |
! .............. ! ! .............. ! .............. !
: ' : ' '
move.w (sin,d0.w*4),d1 * Fetch sine-value.
move.w 2(sin,d0.w*4),d2 * Fetch cosine-vlueOR:
movem.w (sin,d0.w*4),d1-d2 * d1.w:sin(x), d2.w:cos(x)That takes care of the points 1 and 4, but what about point 2, the sine- generator. Certainly Seboz's code couldn't be that easy to beat?
I could have based my routine on Seboz's routine, but then it would have gotten bigger, because it was already size-optimised and with the new sine/cosine- table some extra instructions would even be required..
Finally: I forgot about Seboz's routine and decided to optimise some of the old sine-generators for 68020 instructions, so they would be a bit smaller.
One particular routine was based on the math:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~ :Math alert 2: ~~~~~~~~~~~~~~~~~~~~~~~
sin(à+á)=sin(à)*cos(á)+sin(á)*cos(à)
cos(à+á)=cos(à)*cos(á)+sin(á)*sin(à)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~As you can see, if you know the cos and sin values of the first angle in the array, you can calculate every other one!! Simply let á be the interval between two consecutive angles and we're there.
Ok, now check out the include-file I've made:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ :listing 3: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*==========================================================================
* Sinewave table generator.
* By EarX/~fUn~, 10-5-1998
* 68020 or higher is required!
*==========================================================================
******** GLOBAL CONSTANTS ********
* >WARNING< for these equ's: when using a new 'sintbllen' you must
* recalculate 'cos1' and 'sin1'!
sintbllen: equ 2048 * MUST BE A EXPONENTIAL VALUE OF 2!
sin1: equ 13176774 * sin(2ã/2048)*2^32
cos1: equ 4294947083 * cos(2ã/2048)*2^32
******** MACROS ********
* Macro that returns the modulo of a given angle.
* INPUT: angle: type: data-register (word) or RAM (word)
Do_SinModulo: MACRO angle
andi.w #sintbllen-1,\1
ENDM
* Macro that returns sine & cosine of a given angle.
* PRECONDITION: INIT_SINETABLE has been called!
* INPUT: a0: address of sine_tbl
* base: type: address-register or address or relative address
* inpreg: type: data-register or address-register (lower word)
* contains: angle (0=0ø, sintbllen=360ø)
* OUTPUT: sinreg: type: data-register (long) or address-register
* contains: sine value (signed: -32768 to 32767)
* cosreg: type: data-register (long) or address-register
* contains: cosine value (signed: -32768 to 32767)
Get_SinCos: MACRO base,inpreg,sinreg,cosreg
movem.w (\1,\2.w*4),\3/\4
ENDM
* Macro that returns sine of a given angle.
* PRECONDITION: INIT_SINETABLE has been called!
* INPUT: a0: address of sine_tbl
* base: type: address-register or address or relative address
* inpreg: type: data-register or address-register (lower word)
* contains: angle (0=0ø, sintbllen=360ø)
* OUTPUT: sinreg: type: data-register (long) or address-register
* contains: sine value (signed: -32768 to 32767)
Get_Sin: MACRO base,inpreg,sinreg
move.w (\1,\2.w*4),\3
ENDM
* Macro that returns cosine of a given angle.
* PRECONDITION: INIT_SINETABLE has been called!
* INPUT: a0: address of sine_tbl
* base: type: address-register or address or relative address
* inpreg: type: data-register or address-register (lower word)
* contains: angle (0=0ø, sintbllen=360ø)
* OUTPUT: cosreg: type: data-register (long) or address-register
* contains: cosine value (signed: -32768 to 32767)
Get_Cos: MACRO base,inpreg,sinreg,cosreg
move.w 2(\1,\2.w*4),\3
ENDM
* Creates the a combined sine and cosine table for quick fetching.
* Macro is exactly 96 bytes in length :-)
* INPUT: a0: address of sine_tbl (must be initialized as all zeroes!)
Init_SineTable: MACRO
moveq #$ffffffff,d0 * /cos(0)=1
lsr.l #1,d0 * \(=$7fffffff)
moveq #0,d1 * sin(0)=0
move.l #sin1,d6
move.w #sintbllen/4-1,d7
.genlop:
swap d0 * Get high-word of cosa
swap d1 * Get high-word of sina
move.w d1,2+(sintbllen)*3(a0) * Copy sina in cos-4th quadrant
move.w d0,sintbllen*1(a0) * Copy cosa in sin-2nd quadrant
sub.w d1,2+(sintbllen)*1(a0) * Copy -sina in cos-2nd quadrant
sub.w d0,sintbllen*3(a0) * Copy -cosa in sin-4th quadrant
sub.w d0,2+(sintbllen)*2(a0) * Copy -cosa in cos-3rd quadrant
sub.w d1,sintbllen*2(a0) * Copy -sina in sin-3rd quadrant
move.w d1,(a0)+ * Save sina (16 bit signed value) in first quadrant
move.w d0,(a0)+ * Save cosa (16 bit signed value) in first quadrant
swap d0 * Change cosa back to fixedpoint
swap d1 * Change sina back to fixedpoint
move.l d1,d4 * / Backup sina
move.l d0,d5 * | and cosa
move.l d1,d2 * | for use in
move.l d0,d3 * \ multiplications.
mulu.l d6,d3:d1 * d3:=sin1*sina
mulu.l #cos1,d2:d0 * d2:=cos1*cosa
mulu.l d6,d1:d5 * d0:=sin1*cosa
mulu.l #cos1,d0:d4 * d1:=cos1*sina
sub.l d3,d2 * d2:=(cos1*cosa)-(sin1*sina)
add.l d0,d1 * sina:=(sin1*cosa)+(cos1*sina)
move.l d2,d0 * cosa:=(cos1*cosa)-(sin1*sina)
dbra d7,.genlop
ENDM
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~If you decide to use it you simply do like this:
- Include these MACRO's in your code.
- Put "Init_SineTable" in your code, to initialize the sinetable.
- Put "Get_Sin a0,d0,d1" to get to sine of d0.w in d1.w
Put "Get_Cos a0,d0,d1" to get to cosine of d0.w in d1.w
Put "Get_SinCos a0,d0,d1,d2" to get to sin/cos of d0.w in d1.w/d2.w
Put "Do_SinModulo d0" to get the modulo of the angle in d0.w
Ok, that's it. Blah! I'm tired. Have FUN with this and do you have any comments or suggestions? Please mail me!! >> 1028587@ibk.fnt.hvu.nl
EarX/FUN
