; 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
