#include <time.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define SAMPLERATE 44100 /* samples per second */
#define MIDDLE_C 261.6255653005986346778499935233 /* based on A-440 */
#define SEMITONE 1.0594630943592952645618252949463417 /* 12th root of 2 */
#define STARTFUZZ (SAMPLERATE/100) /* up to 10ms of start time fuzz */
#define EARDIST 127 /* mm from head middle to ear */
#define MACH1 343200 /* mm per second */

#define MAXDELAY ((int)((EARDIST * 2 * SAMPLERATE) / MACH1))

typedef struct thread THREAD;
typedef struct noteparams NOTEPARAMS;
typedef struct playing PLAYING;
typedef struct waveform WAVEFORM;

struct waveform {
  int relprob;
  void *(*init)(PLAYING *);
  double (*sample)(PLAYING *, double);
  void (*done)(void *);
  } ;
#define WAVEFORM_INIT(prob,tag) {	\
	prob,				\
	&wave_init_##tag,		\
	&wave_sample_##tag,		\
	&wave_done_##tag		\
	}

struct noteparams {
  int note;
  int octave;
  int attack;
  int decay;
  const WAVEFORM *wave;
  double loc_r;
  double loc_phi;
  } ;

struct thread {
  unsigned int beatsshift;
  unsigned int beatsamples;
  unsigned int beats;
  NOTEPARAMS **pattern;
  NOTEPARAMS defparams;
  unsigned int beatat;
  unsigned int nextbeat;
  } ;

struct playing {
  int age;
  NOTEPARAMS params;
  void *wavestate;
  double pitch;
  double vol;
  int end;
  int eardelay; /* in samples; +ve => left of centre */
  } ;

static const int notepowers[5] = { 1, 3, 6, 8, 10 };

static int nthreads;
static int athreads;
static THREAD *threads;
static int nplaying;
static int aplaying;
static PLAYING *playing;
static unsigned int age;
static int mburst;

static void init_playing(void)
{
 nplaying = 0;
 aplaying = 8;
 playing = malloc(aplaying*sizeof(*playing));
}

static void init_threads(void)
{
 nthreads = 0;
 athreads = 0;
 threads = 0;
}

static void *wave_init_sine(PLAYING *play __attribute__((__unused__)))
{
 return(0);
}

static double wave_sample_sine(PLAYING *play, double t)
{
 return(1.25*sin((t*2*M_PI*play->pitch)/SAMPLERATE));
}

static void wave_done_sine(void *pv __attribute__((__unused__)))
{
}

static const WAVEFORM waveform_sine = WAVEFORM_INIT(200,sine);

static void *wave_init_square(PLAYING *play __attribute__((__unused__)))
{
 return(0);
}

static double wave_sample_square(PLAYING *play, double t)
{
 return((1&(int)((t*2*play->pitch)/SAMPLERATE))?.3:-.3);
}

static void wave_done_square(void *pv __attribute__((__unused__)))
{
}

static const WAVEFORM waveform_square = WAVEFORM_INIT(100,square);

static void *wave_init_triangle(PLAYING *play __attribute__((__unused__)))
{
 return(0);
}

static double wave_sample_triangle(PLAYING *play, double t)
{
 double v;

 v = (t * play->pitch) / SAMPLERATE;
 v -= (int)v;
 return( (v < .5) ? ((v * 4) - 1) : (3 - (v * 4)) );
}

static void wave_done_triangle(void *pv __attribute__((__unused__)))
{
}

static const WAVEFORM waveform_triangle = WAVEFORM_INIT(100,triangle);

typedef struct priv_rectangle PRIV_RECTANGLE;

struct priv_rectangle {
  double pwm_offset;
  } ;

static void *wave_init_rectangle(PLAYING *play __attribute__((__unused__)))
{
 PRIV_RECTANGLE *r;

 r = malloc(sizeof(PRIV_RECTANGLE));
 r->pwm_offset = (random() & 65535) / 65536.0;
 return(r);
}

static double wave_sample_rectangle(PLAYING *play, double t)
{
 PRIV_RECTANGLE *r;
 double pv;
 double ev;

 r = play->wavestate;
 pv = (t * play->pitch) / SAMPLERATE;
 pv -= (int)pv;
 ev = ((t * 2) / SAMPLERATE) + r->pwm_offset;
 ev -= (int)ev;
 if (ev > .5) ev = 1 - ev;
 return((pv<ev)?.3:-.3);
}

static void wave_done_rectangle(void *pv)
{
 free(pv);
}

static const WAVEFORM waveform_rectangle = WAVEFORM_INIT(50,rectangle);

static const WAVEFORM *waveforms[]
 = { &waveform_sine,
     &waveform_square,
     &waveform_triangle,
     &waveform_rectangle,
     0 };

static void start_note(NOTEPARAMS *def, NOTEPARAMS *np)
{
 PLAYING *p;
 double locx;
 double locy;

 if (nplaying >= aplaying) playing = realloc(playing,(aplaying=nplaying+8)*sizeof(*playing));
 p = &playing[nplaying];
 p->age = - (random() % STARTFUZZ) - MAXDELAY;
 p->params = *def;
 if (np->note >= 0) p->params.note = np->note;
 if (np->octave >= 0) p->params.octave = np->octave;
 if (np->attack >= 0) p->params.attack = np->attack;
 if (np->decay >= 0) p->params.decay = np->decay;
 if (np->wave) p->params.wave = np->wave;
 if (np->loc_r >= 0) p->params.loc_r = np->loc_r;
 if (np->loc_phi >= 0) p->params.loc_phi = np->loc_phi;
 p->pitch = MIDDLE_C * pow(SEMITONE,(np->octave*12)+notepowers[np->note]);
 p->wavestate = (*p->params.wave->init)(p);
 p->vol = log(MIDDLE_C) / log(p->pitch);
 p->vol *= p->vol;
 p->end = p->params.attack + p->params.decay + 2;
 locx = p->params.loc_r * cos(p->params.loc_phi);
 locy = p->params.loc_r * sin(p->params.loc_phi);
 p->eardelay = (hypot(locx-EARDIST,locy) - hypot(locx+EARDIST,locy)) * (SAMPLERATE / MACH1);
 nplaying ++;
}

static void start_notes(void)
{
 int i;
 THREAD *t;
 NOTEPARAMS *p;

 for (i=nthreads-1;i>=0;i--)
  { t = &threads[i];
    if (age >= t->beatat)
     { p = t->pattern[t->nextbeat];
       if (p) start_note(&t->defparams,p);
       t->beatat += t->beatsamples;
       t->nextbeat ++;
       if (t->nextbeat >= t->beats) t->nextbeat = 0;
     }
  }
}

static void gen_sample(void)
{
 double ls;
 double rs;
 int is;
 int i;
 PLAYING *p;
 int live;

 double notesample(double t)
  { double v;
    if (t < 0) return(0);
    if (t < p->params.attack)
     { v = 1 - cos((t*M_PI)/p->params.attack);
     }
    else if (t < p->end)
     { v = 1 + cos(((t-p->params.attack)*M_PI)/p->params.decay);
     }
    else
     { live --;
       return(0);
     }
    return(v*p->vol*(*p->params.wave->sample)(p,t));
  }

 ls = 0;
 rs = 0;
 for (i=nplaying-1;i>=0;i--)
  { p = &playing[i];
    live = 2;
    ls += notesample(p->age);
    rs += notesample(p->age-p->eardelay);
    if (live < 1)
     { nplaying --;
       (*p->params.wave->done)(p->wavestate);
       if (i != nplaying) playing[i] = playing[nplaying];
       continue;
     }
    p->age ++;
  }
 is = (ls * 2000) + 32768;
 if (is < 0) is = 0; else if (is > 65535) is = 65535;
 putchar(is&0xff);
 putchar((is>>8)^0x80);
 is = (rs * 2000) + 32768;
 if (is < 0) is = 0; else if (is > 65535) is = 65535;
 putchar(is&0xff);
 putchar((is>>8)^0x80);
}

static const WAVEFORM *pick_waveform(void)
{
 int p;
 int i;

 p = 0;
 for (i=0;waveforms[i];i++) p += waveforms[i]->relprob;
 p = random() % p;
 for (i=0;;i++)
  { p -= waveforms[i]->relprob;
    if (p < 0) return(waveforms[i]);
  }
 abort();
}

static void pick_envelope(NOTEPARAMS *p)
{
 switch (random() & 7)
  { case 0: case 1: case 2:
       p->attack = SAMPLERATE / 100;
       p->decay = SAMPLERATE / 4;
       break;
    case 3:
       p->attack = SAMPLERATE / 4;
       p->decay = SAMPLERATE / 100;
       break;
    case 4:
       p->attack = SAMPLERATE / 100;
       p->decay = SAMPLERATE / 2;
       break;
    case 5:
       p->attack = SAMPLERATE / 2;
       p->decay = SAMPLERATE / 100;
       break;
    case 6:
       p->attack = SAMPLERATE / 50;
       p->decay = SAMPLERATE;
       break;
    case 7:
       p->attack = SAMPLERATE / 5;
       p->decay = SAMPLERATE;
       break;
  }
}

static void mutate_start_or_kill(void)
{
 int n;
 int i;
 THREAD *t;
 double x;
 double y;

 n = random() & 15;
 if (n > nthreads)
  { if (nthreads >= athreads) threads = realloc(threads,(athreads=nthreads+8)*sizeof(*threads));
    t = &threads[nthreads];
    t->beatsshift = 12 + (random() & 3);
    t->beatsamples = 1 << t->beatsshift;
    t->beats = (2 + (random() & 3)) << (random() & 3);
    t->pattern = malloc(t->beats*sizeof(*t->pattern));
    for (i=t->beats-1;i>=0;i--) t->pattern[i] = 0;
    t->defparams.note = random() % 5;
    switch (random() % 10)
     { case 0:
	  t->defparams.octave = -2;
	  break;
       case 1: case 2:
	  t->defparams.octave = -1;
	  break;
       case 3: case 4: case 5: case 6:
	  t->defparams.octave = 0;
	  break;
       case 7: case 8:
	  t->defparams.octave = 1;
	  break;
       case 9:
	  t->defparams.octave = 2;
	  break;
     }
    pick_envelope(&t->defparams);
    t->defparams.wave = pick_waveform();
    x = ((random() & 65535) - 32768.5) / 16;
    y = ((random() & 65535) / 16.0) + EARDIST;
    t->defparams.loc_r = hypot(x,y);
    t->defparams.loc_phi = atan2(y,x);
    t->pattern[0] = malloc(sizeof(*t->pattern[0]));
    *t->pattern[0] = t->defparams;
    t->beatat = ((age / t->beatsamples) + 1) * t->beatsamples;
    t->nextbeat = 0;
    nthreads ++;
  }
 else if (n < nthreads)
  { nthreads --;
    t = &threads[n];
    for (i=t->beats-1;i>=0;i--) free(t->pattern[i]);
    free(t->pattern);
    if (n < nthreads) threads[n] = threads[nthreads];
  }
}

static void mutate_voice_or_silent(void)
{
 int silentbeats;
 int voicedbeats;
 int tx;
 THREAD *t;
 int px;
 int n;

 silentbeats = 0;
 voicedbeats = 0;
 for (tx=nthreads-1;tx>=0;tx--)
  { t = &threads[tx];
    for (px=t->beats-1;px>=0;px--)
     { if (t->pattern[px]) voicedbeats ++; else silentbeats ++;
     }
  }
 n = random() % (silentbeats + voicedbeats + 1);
 if (n < voicedbeats)
  { for (tx=nthreads-1;tx>=0;tx--)
     { t = &threads[tx];
       for (px=t->beats-1;px>=0;px--)
	{ if (t->pattern[px])
	   { if (n < 1)
	      { free(t->pattern[px]);
		t->pattern[px] = 0;
		return;
	      }
	     n --;
	   }
	}
     }
  }
 else if (n > voicedbeats)
  { n -= voicedbeats;
    for (tx=nthreads-1;tx>=0;tx--)
     { t = &threads[tx];
       for (px=t->beats-1;px>=0;px--)
	{ if (! t->pattern[px])
	   { if (n < 1)
	      { t->pattern[px] = malloc(sizeof(*t->pattern[px]));
		*t->pattern[px] = t->defparams;
		return;
	      }
	     n --;
	   }
	}
     }
  }
}

static void mutate_half_beats(void)
{
 THREAD *t;
 NOTEPARAMS **pv;
 int i;

 t = &threads[random()%nthreads];
 if (random() & 1)
  { if (t->beatsshift > 11)
     { pv = malloc(t->beats*2*sizeof(*pv));
       for (i=t->beats-1;i>=0;i--)
	{ pv[i+i+1] = 0;
	  pv[i+i] = t->pattern[i];
	}
       free(t->pattern);
       t->pattern = pv;
       t->beats *= 2;
       t->beatsshift --;
       t->beatsamples >>= 1;
       t->nextbeat <<= 1;
     }
  }
 else
  { if ((t->beatsshift < 15) & !(t->beats & 1))
     { pv = malloc((t->beats>>1)*sizeof(*pv));
       for (i=(t->beats>>1)-1;i>=0;i--)
	{ pv[i] = t->pattern[i+i];
	  free(t->pattern[i+i+1]);
	}
       free(t->pattern);
       t->pattern = pv;
       t->beats >>= 1;
       t->beatsshift ++;
       t->beatsamples <<= 1;
       t->nextbeat >>= 1;
     }
  }
}

static void mutate_parameters(NOTEPARAMS *p)
{
 int n;
 int r;

 n = p->note + (5 * p->octave);
 r = (random() % 25) - 10;
 if (r < n) n --; else if (r > n) n ++;
 p->octave = n / 5;
 p->note = n % 5;
 if (random() & 1) pick_envelope(p);
 if (random() & 1) p->wave = pick_waveform();
}

static void mutate_change_note(void)
{
 int tx;
 THREAD *t;
 int px;
 int c;
 NOTEPARAMS *found;

 c = 0;
 for (tx=nthreads-1;tx>=0;tx--)
  { t = &threads[tx];
    for (px=t->beats-1;px>=0;px--)
     { if (t->pattern[px])
	{ c ++;
	  if (! (random() % c))
	   { found = t->pattern[px];
	   }
	}
     }
  }
 if (! c) return;
 mutate_parameters(found);
}

static void mutate_change_default(void)
{
 mutate_parameters(&threads[random()%nthreads].defparams);
}

/*
 * Mutations possible:
 *
 * - Start new thread
 * - Kill existing thread
 * - Change an existing note parameter
 * - Change an existing thread default parameter
 * - Silence a thread beat
 * - Voice a silent thread beat
 * - Divide/multiply thread pattern length by 2/3
 * - Introduce half beats
 * - Drop half beats
 *
 * Probabilities change depending on existing state.
 */
static void mutate(void)
{
 if (! nthreads)
  { mutate_start_or_kill();
    return;
  }
 switch (random() & 63)
  { case 0:
       mutate_start_or_kill();
       break;
    case 1: case 2: case 3: case 4:
       mutate_voice_or_silent();
       break;
    case 5:
       mutate_half_beats();
       break;
    case 6: case 7: case 8:
       mutate_change_note();
       break;
    case 9:
       mutate_change_default();
       break;
    case 63:
       while (random() & 15) mburst ++;
       break;
  }
}

int main(void);
int main(void)
{
 srandom(time(0));
 init_playing();
 init_threads();
 mburst = 44100;
 age = 0;
 while (1)
  { start_notes();
    gen_sample();
    if (mburst > 0)
     { mutate();
       mburst --;
     }
    else if (! (random() & 32767))
     { mutate();
     }
    age ++;
  }
}