#include <stdio.h>
#include <stdint.h>
#include <math.h>
#include <string.h>

#define RATE 16000.0

float waveform[2048];

typedef struct oscillator {
	uint32_t phasor;
	uint32_t step;
	float value;
} osc_t;

void update_osc(osc_t* o) {
	o->phasor += o->step;
	o->value = waveform[o->phasor >> (32-11)];
}

void set_osc_freq(osc_t* o, float f) {
	o->step = (uint32_t) ((f / RATE) * UINT32_MAX);
}

int main() {
	osc_t oscbank[16];
	float out, x;
	int16_t intout;

	// Sine wave waveshaped with T2(x) + T3(x)
	// https://en.wikipedia.org/wiki/Chebyshev_polynomials
	for(int i=0; i<2048; i++) {
		x = (float) sin(i*2*M_PI/2048);
		waveform[i] = (4*x*x*x + 2*x*x -3*x - 1) / 2;
	}

	// Four detuned A major (augmented) chords
	set_osc_freq(&oscbank[0], 110);
	set_osc_freq(&oscbank[1], 110 + .01);
	set_osc_freq(&oscbank[2], 110 + .025);
	set_osc_freq(&oscbank[3], 110 + .05);

	set_osc_freq(&oscbank[4], 138.59);
	set_osc_freq(&oscbank[5], 138.59 + .025);
	set_osc_freq(&oscbank[6], 138.59 + .050);
	set_osc_freq(&oscbank[7], 138.59 + .075);

	set_osc_freq(&oscbank[8],  164.81);
	set_osc_freq(&oscbank[9],  164.81 + .075);
	set_osc_freq(&oscbank[10], 164.81 + .1);
	set_osc_freq(&oscbank[11], 164.81 + .2);

	set_osc_freq(&oscbank[12], 207.65 + 0.1);
	set_osc_freq(&oscbank[13], 207.65 + 0.15);
	set_osc_freq(&oscbank[14], 207.65 + 0.2);
	set_osc_freq(&oscbank[15], 207.65 + 0.3);

	while(1) {
		out = 0;
		for(int i=0; i<16; i++) {
			update_osc(&oscbank[i]);
			out += (1.0/16.0)*oscbank[i].value;
		}
		intout = (int16_t) (INT16_MAX)*out;
		fwrite(&intout, 2, 1, stdout);
	}

	return 1;
}