#include <u.h>
#include <libc.h>
typedef struct Complex Complex;
struct Complex{double re, im;};
double sine[512];
Complex
mul(Complex a, Complex b){
	double re;
	re = a.re * b.re - a.im * b.im;
	a.im = a.re * b.im + a.im * b.re;
	a.re = re;
	return a;
}
Complex
div(Complex a, Complex b){
	double re, mag;
	mag = b.re * b.re + b.im * b.im;
	re = (a.re * b.re + a.im * b.im) / mag;
	a.im = (a.im * b.re - a.re * b.im) / mag;
	a.re = re;
	return a;
}
Complex
oneminus(Complex a){
	a.re = 1.0 - a.re;
	a.im = 1.0 - a.im;
	return a;
}
void
main(int, char**){
	uint uφ, vφ, uφincr, vφincr, i;
	Complex u, v, x;
	double w = 0.9, norm;
	uchar buf[44100 * 5 * 4];
	int fd, n, N = 9, s;
	fd = open("/dev/audio", OWRITE);
	if(fd < 0) sysfatal("open: %r");
	for(n = 0; n < 512; n++)
		sine[n] = sin(n * PI / 256.0);
	norm = (1.0 - pow(w, N + 1)) / (1.0 - w);
	uφ = vφ = 0;
	uφincr = (512 * 440 << 12) / 44100;
	vφincr = (512 * 110 << 12) / 44100;
	for(n = 0; n < sizeof buf; n += 4){
		i = uφ >> 12;
		uφ += uφincr;
		u.re = sine[i + 256 & 511];
		u.im = sine[i & 511];
		i = vφ >> 12;
		vφ += vφincr;
		x.re = v.re = w * sine[i + 256 & 511];
		x.im = v.im = w * sine[i & 511];
		for(i = 0; i < N; i++)
			x = mul(x, v);
		x = div(mul(u, oneminus(x)), oneminus(v));
		s = (0.7 * x.im + 0.3 * x.re) / norm * 0x2000;
		buf[n] = s;
		buf[n | 1] = s >> 8;
		s = (0.7 * x.re + 0.3 * x.im) / norm * 0x2000;
		buf[n | 2] = s;
		buf[n | 3] = s >> 8;
	}
	write(fd, buf, sizeof buf);
	exits(nil);
}