; FILE: Source:sinegen.ASM REV: 31 --- 16-bit sinetable generator ; History ; 31 18th September 1998: 1st version. ; IFGT 0 Inspiration for this document and source came from PAC/#amycoders who needed good&short sinetable generator. My friend ArtDent coded this kind of routine years ago, but unfortunately he didn't backup his amiga sources when he went pc. Anyways he still remembered the principle well and he pointed me the algorithm to use. This whole document and source was written by me (Piru) in 5 hours. sine&cosine table generation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Lets have a look at sine and cosine graph: pi 2pi _ | | |/|\| | | --/-+-\-+-/-- | | |\|/| 0 | T | | 1/2pi 3/2pi pi 3/2pi _ | | _ |\| | |/| --+-\-+-/-+-- | |\_/| | 0 | | 1/2pi 2pi We notice that sine is phase shifted 90 degrees compared to cosine. Also we notice that both sine and cosine are symmetrical to 1/2pi and pi, thus can be easily mirrored. So we need to calculate only 90 degrees of either sine or cosine and we can derive whole table from it and also the other function. These are the formulas to calculate sin x and cos x: sin x = x - x^3 / 3! + x^5 / 5! - x^7 / 7! + ... cos x = 1 - x^2 / 2! + x^4 / 4! - x^6 / 6! + ... x is real, 0 <= x <= 1/2pi Out of these two the latter (cos x) is easier to calculate. You can save space by combining sine and cosine tables. Just take last 90 degrees of cosine before cosine table and you have sinetable at table - 90 degrees. :) So after thinking a while I came up with this pseudocode routine that calculates 90 degrees of sine + 360 degrees cosine: in: table, tablesize (90 degrees * 5) quart = tablesize / 5 x = 0; x_add = (1/2 * pi) / quart for q = 0 to (quart - 1) fact = 1; d = 0; cosx = 1; powx = 1 powx_mul = - (x * x) ; rem this will magically toggle sign repeat powx = powx * powx_mul d++; fact = fact * d d++; fact = fact * d cosx = cosx + powx / fact until d = 12 table[quart - q] = cosx ; rem /¯ table[quart + q] = cosx ; rem ¯\ table[quart * 3 - q] = -cosx ; rem \_ table[quart * 3 + q] = -cosx ; rem _/ table[quart * 5 - q] = cosx ; rem /¯ x = x + x_add endfor Then I just coded this in 020+ asm adding fixedpoint math and stuff: ENDC TESTSINE SET 0 IFNE TESTSINE Main lea (sine,pc),a0 move.l #256,d0 bsr sinegen rts sine ds.w 256 cosine ds.w 256*4 ENDC ; 68020+ 16:16 fixedpoint sinetable generator. ; Coded by Harry "Piru" Sintonen. ; Not specially optimized as usually this thing is ran only once at ; init time. 68060 will woe on 64 bit muls & swaps - who cares ;) ; IN: a0.l=pointer to array of word (will contain 450 degree 16-bit sinetable) ; d0.l=wordsper90degrees ; OUT: d0.l=0 sinegen movem.l d0-d7/a0-a5,-(sp) move.l #26353589,d1 ; pi/2*65536*256 divu.l d0,d1 move.l d1,a5 add.l d0,d0 add.l d0,a0 lea 0(a0,d0.l*2),a2 lea 0(a0,d0.l*4),a4 move.l a0,a1 move.l a2,a3 addq.l #2,a1 ; these two can be removed addq.l #2,a2 ; really ;) moveq #0,d0 ; x moveq #12,d7 .oloop move.l d0,d5 moveq #1,d1 lsr.l #8,d5 swap d1 ; 1<<16 = cos x move.l d1,d3 mulu.l d5,d4:d5 move.w d4,d5 moveq #0,d2 ; d swap d5 moveq #1,d6 ; factorial neg.l d5 ; change sign of powx .iloop muls.l d5,d4:d3 ; calculate x^d move.w d4,d3 swap d3 move.l d3,d4 addq.l #1,d2 ; calculate d! mulu.l d2,d6 addq.l #1,d2 mulu.l d2,d6 divs.l d6,d4 add.l d4,d1 ; cos x += x^d / d! cmp.l d7,d2 bne.b .iloop lsr.l #1,d1 tst.w d1 ; if d1=$8000 then d1=d1-1 ;) dbpl d1,.rule .rule move.w d1,(a0)+ move.w d1,-(a1) move.w d1,-(a4) neg.w d1 move.w d1,-(a2) move.w d1,(a3)+ add.l a5,d0 subq.l #1,(sp) ; watch out - don't mess with stack:) bne.b .oloop movem.l (sp)+,d0-d7/a0-a5 rts