sin(x - y) + sin(x + y) = 2 * sin x * cos y <=>
sin(x + y) = 2 * sin x * cos y - sin(x - y)

k = 2 * cos d

sin[i + d] = sin[i] * k - sin[i - d]

table	= $0400		;; preferably a low page

loop	lda #<index
	ldy #>index
	jsr $bba2
	
	jsr $e277

	lda $61
	adc #7
	sta $61

	jsr $bc9b
	lda $65

index	.byte $90
	.byte $00
	sta table

	inc *-2
	bne loop

table	= $0400		;; Preferably a low page. Must be paged aligned!

loop	lda #<index	;; Load 5-byte float at 'index' into FAC, the fraction
	ldy #>index	;; of which is stepped between 0/256..255/256.
	jsr $bba2	;; However an integer bias is also added in order to fix
			;; the exponent and make hence it possible to increment the
			;; fraction as a normal binary byte, e.g. a version of the
			;; classic x86 float-to-int conversion trick.
	
	jsr $e277	;; Now calculate sine of FAC. Except skip the initial part
			;; of the BASIC function which divides by 2*PI to get
			;; a fraction out of radians since we've already got one.
			;; The integer bias is taken care by BASIC since sin()
			;; is supposed to be periodic.

	lda #<bias	;; Convert the output in FAC from a float in -1..+1 to a
	ldy #>bias	;; fixed-point value in -128..+127 at the LSB of the
	jsr $b867	;; mantissa by employing the same trick as before of
	lda $65		;; adding a high integer bias.

index	.byte $90	;; This is both a float *and* a piece of code. The exponent
	.byte $00	;; ($90 corresponds to 2^16) fixes our 8-bit fraction as the
	sta table	;; third byte of the mantissa and STA address' LSB (don't
			;; forget that BASIC floats are big-endian!). And $90/$00
			;; interpreted as code simply correspond to a harmless BCC *+2.
			;; Note that while the STA's opcode is an integer part which
			;; shouldn't affect the result, the table's high byte *does*
			;; serve as a small offset shifting the result by up to one
			;; index value. Some might even argue that placing the table
			;; at $8000 would produce the 'proper' rounding.

	inc *-2
	bne loop

bias	= $befa		;; A float with an exponent of $99 (2^25) and an LSB of
;;	.byte $99	;; zero is used to convert the output to binary. Such byte
;;	.byte $02	;; sequences can be found in six places in the BASIC/Kernal
;;	.byte $01	;; ROMs, at $befa/$bf04/$bf09/$fd53/$fd56/$ff38.
;;	.byte $a9	;; A version with an LSB of $80 would have been useful to
;;	.byte $00	;; create unsigned output (e.g. between $00 and $ff with the
			;; origin at $80) but unfortunately doesn't seem to exist.

0 fori=0to255:poke1024+i,sin(i*.02454)*128and255:next

	ldx#127
prefill
	sta cosine,x
	lda#19
	sta cosine+52,x
	lda#250
	sta cosine+128+52,x
	dex
	bpl prefill

	ldy#30
	clc
smoothloop
	lda cosine,x
	adc cosine,y
	ror
	sta cosine,x
	iny
	inx
	bne smoothloop
	dey
	bne smoothloop

import numpy as N

def sim(si,jo,dco,lv,hv,pw):
	v=N.zeros((256,),N.int)+lv
	v[dco:][:pw]*=0
	v[dco:][:pw]+=hv

	nv=1

	i=si
	j=(i+jo)&255
	while 1:
		nv=(v[i]+v[j])+(nv&1)
		v[i]=nv/2
		i=(i+1)&255
		j=(j+1)&255
		if(i&255==0):
			j=(j-1)&255
			if j==0:
				break

	v=(v+N.take(v,(N.arange(256)+64)%256))*256
	v-=N.sum(v)/256
	v=N.cumsum(v)/256
	v-=N.sum(v)/256
	v[0]+=0xfffe00
	v=N.cumsum(v)/256
	return v

v=sim(255, 31, 81, 0, 13097, 131)