#include "amy.h"


// For checking assumptions about bitwidths.
#include <assert.h>

#include "sine_lutset_fxpt.h"
#include "saw_lutset_fxpt.h"
#include "triangle_lutset_fxpt.h"

// For hardware random on ESP
#ifdef ESP_PLATFORM
#include <esp_system.h>
#endif

// Multiply used to apply amplitude to unit waveforms. sample is first arg (<=1), amp is 2nd.
//#define MULA_SS MUL6A_SS
#define MULA_SS(a, b) MUL5A_SS(a, b)


/* Dan Ellis libblosca functions */
const LUT *choose_from_lutset(float period, const LUT *lutset) {
    // Select the best entry from a lutset for a given period. 
    //
    // Args:
    //    period: (float) Target period of waveform in fractional samples.
    //    lutset: Sorted list of LUTs, as generated by create_lutset().
    //
    // Returns:
    //   One of the LUTs from the lutset, best suited to interpolating to generate
    //   a waveform of the desired period.
    // Use the earliest (i.e., longest, most harmonics) LUT that works
    // (i.e., will not actually cause aliasing).
    // So start with the highest-bandwidth (and longest) LUTs, but skip them
    // if they result in aliasing.
    const LUT *lut_table = NULL;
    int lut_size = 0;
    int lut_index = 0;
    while(lutset[lut_index].table_size > 0) {
        lut_table = &lutset[lut_index];
        lut_size = lutset[lut_index].table_size;
        // What proportion of nyquist does the highest harmonic in this table occupy?
        float lut_bandwidth = 2 * lutset[lut_index].highest_harmonic / (float)lut_size;
        // To complete one cycle of <lut_size> points in <period> steps, each step
        // will need to be this many samples:
        float lut_hop = lut_size / period;
        // If we have a signal with a given bandwidth, but then speed it up by 
        // skipping lut_hop samples per sample, its bandwidth will increase 
        // proportionately.
        float interp_bandwidth = lut_bandwidth * lut_hop;
        //printf("period=%f freq=%f lut_size=%d interp_bandwidth=%f\n", period, ((float)AMY_SAMPLE_RATE)/period, lut_size, interp_bandwidth);
        if (interp_bandwidth < 0.9) {
            // No aliasing, even with a 10% buffer (i.e., 19.8 kHz).
            break;
        }
        ++lut_index;
    }
    // At this point, we either got to the end of the LUT table, or we found a
    // table we could interpolate without aliasing.

    return lut_table;
}

// Multiple versions of render_lut with different features to avoid branches in sample loop.

#define RENDER_LUT_PREAMBLE \
    if(lut == NULL)  return phase;\
    int lut_mask = lut->table_size - 1; \
    int lut_bits = lut->log_2_table_size; \
    SAMPLE sample = 0; \
    SAMPLE max_value = 0; \
    SAMPLE current_amp = incoming_amp;                                      \
    SAMPLE incremental_amp = SHIFTR(ending_amp - incoming_amp, BLOCK_SIZE_BITS);

#define MOD_PART_MOD  \
            total_phase += S2P(mod[i]);

// Feedback is taken before output scaling.
#define FEEDBACK_PART_FB  \
            past1 = past0;  \
            past0 = sample;  \
            total_phase += S2P(MUL4_SS(feedback_level, SHIFTR(past1 + past0, 1)));

#define RENDER_LUT_GUTS(MOD_PART, FEEDBACK_PART, INTERP_PART) \
            MOD_PART \
            FEEDBACK_PART \
            int16_t base_index = INT_OF_P(total_phase, lut_bits); \
            SAMPLE frac = S_FRAC_OF_P(total_phase, lut_bits); \
            SAMPLE b = L2S(lut->table[base_index]); \
            SAMPLE c = L2S(lut->table[(base_index + 1) & lut_mask]); \
            INTERP_PART

#define INTERP_LINEAR \
            sample = b + MUL0_SS(c - b, frac);

// Miller's optimization -
// https://github.com/pure-data/pure-data/blob/db777311d808bb3ba728b94ab067f8d333b7d0c2/src/d_array.c#L831C1-L833C76
// outlet_float(x->x_obj.ob_outlet, b + frac * (
//    cminusb - 0.1666667f * (1.-frac) * (
//        (d - a - 3.0f * cminusb) * frac + (d + 2.0f*a - 3.0f*b))));#
#define INTERP_CUBIC \
            SAMPLE a = L2S(lut->table[(base_index - 1) & lut_mask]); \
            SAMPLE d = L2S(lut->table[(base_index + 2) & lut_mask]); \
            SAMPLE cminusb = c - b; \
            SAMPLE fr_d_ma_m3cmb = MUL0_SS(d - a - cminusb - SHIFTL(cminusb, 1), frac); \
            SAMPLE next_bit = MUL0_SS(fr_d_ma_m3cmb + d + SHIFTL(a - b, 1) - b, MUL0_SS(F2S(1.0f) - frac, F2S(0.16666666666667f))); \
            sample = b + MUL0_SS(cminusb - next_bit, frac);

#define RENDER_LUT_LOOP_END \
            SAMPLE value = buf[i] + MULA_SS(sample, current_amp);	\
            buf[i] = value;                            \
            if (value < 0) value = -value;             \
            if (value > max_value) max_value = value;  \
            current_amp += incremental_amp; \
            phase = P_WRAPPED_SUM(phase, step);


#define NOTHING ;


PHASOR render_lut_fm_fb(SAMPLE* buf,
                        PHASOR phase, 
                        PHASOR step,
                        SAMPLE incoming_amp, SAMPLE ending_amp,
                        const LUT* lut,
                        SAMPLE* mod, SAMPLE feedback_level, SAMPLE* last_two, SAMPLE *pmax_value) { 
    AMY_PROFILE_START(RENDER_LUT_FM_FB)
    RENDER_LUT_PREAMBLE
    SAMPLE past0 = 0, past1 = 0;
    sample = last_two[0];
    past0 = last_two[1];
    for(uint16_t i = 0; i < AMY_BLOCK_SIZE; i++) {
        PHASOR total_phase = phase;

        RENDER_LUT_GUTS(MOD_PART_MOD, FEEDBACK_PART_FB, INTERP_LINEAR)

        RENDER_LUT_LOOP_END
    }
    last_two[0] = sample;
    last_two[1] = past0;
    *pmax_value = max_value;
    AMY_PROFILE_STOP(RENDER_LUT_FM_FB)
    return phase;
}

PHASOR render_lut_fb(SAMPLE* buf,
                     PHASOR phase,
                     PHASOR step,
                     SAMPLE incoming_amp, SAMPLE ending_amp,
                     const LUT* lut,
                     SAMPLE feedback_level, SAMPLE* last_two, SAMPLE *pmax_value) {
    AMY_PROFILE_START(RENDER_LUT_FB)
    RENDER_LUT_PREAMBLE
    SAMPLE past0 = 0, past1 = 0;
    sample = last_two[0];
    past0 = last_two[1];
    for(uint16_t i = 0; i < AMY_BLOCK_SIZE; i++) {
        PHASOR total_phase = phase;

        RENDER_LUT_GUTS(NOTHING, FEEDBACK_PART_FB, INTERP_LINEAR)

        RENDER_LUT_LOOP_END
    }
    last_two[0] = sample;
    last_two[1] = past0;
    *pmax_value = max_value;
    AMY_PROFILE_STOP(RENDER_LUT_FB)
    return phase;
}

PHASOR render_lut_fm(SAMPLE* buf,
                     PHASOR phase,
                     PHASOR step,
                     SAMPLE incoming_amp, SAMPLE ending_amp,
                     const LUT* lut,
                     SAMPLE* mod,
                     SAMPLE* pmax_value) {
    AMY_PROFILE_START(RENDER_LUT_FM)
    RENDER_LUT_PREAMBLE
    for(uint16_t i = 0; i < AMY_BLOCK_SIZE; i++) {
        PHASOR total_phase = phase;

        RENDER_LUT_GUTS(MOD_PART_MOD, NOTHING, INTERP_LINEAR)

        RENDER_LUT_LOOP_END
    }
    *pmax_value = max_value;
    AMY_PROFILE_STOP(RENDER_LUT_FM)
    return phase;
}

PHASOR render_lut(SAMPLE* buf,
                  PHASOR phase,
                  PHASOR step,
                  SAMPLE incoming_amp, SAMPLE ending_amp,
                  const LUT* lut,
                  SAMPLE* pmax_value) {
    AMY_PROFILE_START(RENDER_LUT)
    RENDER_LUT_PREAMBLE
    for(uint16_t i = 0; i < AMY_BLOCK_SIZE; i++) {
        PHASOR total_phase = phase;

        RENDER_LUT_GUTS(NOTHING, NOTHING, INTERP_LINEAR)

        RENDER_LUT_LOOP_END
     }
    *pmax_value = max_value;
    AMY_PROFILE_STOP(RENDER_LUT)
    return phase;
}

PHASOR render_lut_cub(SAMPLE* buf,
                      PHASOR phase,
                      PHASOR step,
                      SAMPLE incoming_amp, SAMPLE ending_amp,
                      const LUT* lut,
                      SAMPLE *pmax_value) {
    AMY_PROFILE_START(RENDER_LUT_CUB)
    RENDER_LUT_PREAMBLE
    for(uint16_t i = 0; i < AMY_BLOCK_SIZE; i++) {
        PHASOR total_phase = phase;

        RENDER_LUT_GUTS(NOTHING, NOTHING, INTERP_CUBIC)

        RENDER_LUT_LOOP_END
    }
    *pmax_value = max_value;
    AMY_PROFILE_STOP(RENDER_LUT_CUB)
    return phase;
}

/* Audio in */

void audio_in_note_on(uint16_t osc, uint8_t channel) {
    // do i need to do anything here? probably not
}

void external_audio_in_note_on(uint16_t osc, uint8_t channel) {
    // do i need to do anything here? probably not
}

SAMPLE render_audio_in(SAMPLE * buf, uint16_t osc, uint8_t channel) {
    uint16_t c = 0;
    for(uint16_t i=channel;i<AMY_BLOCK_SIZE*AMY_NCHANS;i=i+(AMY_NCHANS)) {
        buf[c++] = SMULR7(L2S(amy_in_block[i]), F2S(msynth[osc]->amp));
    }
    // We have to return something for max_value or else the zero-amp reaper will come along. 
    return F2S(1.0); //max_value;
}

SAMPLE render_external_audio_in(SAMPLE *buf, uint16_t osc, uint8_t channel) {
    uint16_t c = 0;
    for(uint16_t i=channel;i<AMY_BLOCK_SIZE*AMY_NCHANS;i=i+(AMY_NCHANS)) {
        buf[c++] = SMULR7(L2S(amy_external_in_block[i]), F2S(msynth[osc]->amp));
    }
    // We have to return something for max_value or else the zero-amp reaper will come along. 
    return F2S(1.0); //max_value;
}


/* Pulse wave */
void pulse_note_on(uint16_t osc, float freq) {
    // Moved to lazy.
    synth[osc]->lut = NULL;
}

void _pulse_note_on(uint16_t osc) {
    //printf("pulse_note_on: time %lld osc %d logfreq %f amp %f last_amp %f\n", amy_global.total_blocks*AMY_BLOCK_SIZE, osc, synth[osc]->logfreq, msynth[osc]->amp, msynth[osc]->last_amp);
    if (synth[osc]->lut == NULL) {
        float freq = freq_of_logfreq(msynth[osc]->logfreq);
        float period_samples = (float)AMY_SAMPLE_RATE / freq;
        synth[osc]->lut = choose_from_lutset(period_samples, saw_fxpt_lutset);
    }
}

SAMPLE render_lpf_lut(SAMPLE* buf, uint16_t osc, int8_t is_square, int8_t direction, SAMPLE dc_offset) {
    AMY_PROFILE_START(RENDER_LPF_LUT)
    // Common function for pulse and saw.
    float freq = freq_of_logfreq(msynth[osc]->logfreq);
    PHASOR step = F2P(freq / (float)AMY_SAMPLE_RATE);  // cycles per sec / samples per sec -> cycles per sample
    SAMPLE amp = direction * F2S(msynth[osc]->amp);
    SAMPLE last_amp = direction * F2S(msynth[osc]->last_amp);
    PHASOR pwm_phase = synth[osc]->phase;
    SAMPLE max_value;
    synth[osc]->phase = render_lut_cub(buf, synth[osc]->phase, step, last_amp, amp, synth[osc]->lut, &max_value);
    if (is_square) {  // For pulse only, add a second delayed negative LUT wave.
        float duty = msynth[osc]->duty;
        if (duty < 0.01f) duty = 0.01f;
        if (duty > 0.99f) duty = 0.99f;
        pwm_phase = P_WRAPPED_SUM(pwm_phase, F2P(msynth[osc]->last_duty));
        // Second pulse is given some blockwise-constant FM to maintain phase continuity across blocks.
        PHASOR delta_phase_per_sample = F2P((duty - msynth[osc]->last_duty) / AMY_BLOCK_SIZE);
        render_lut_cub(buf, pwm_phase, step + delta_phase_per_sample, -last_amp, -amp, synth[osc]->lut, &max_value);
        msynth[osc]->last_duty = duty;
    }
    // Remember last_amp.
    msynth[osc]->last_amp = msynth[osc]->amp;
    AMY_PROFILE_STOP(RENDER_LPF_LUT)
    return max_value;
}

SAMPLE render_pulse(SAMPLE* buf, uint16_t osc) {
    // Second (negative) impulse is <duty> cycles later.
    _pulse_note_on(osc);
    return render_lpf_lut(buf, osc, true, 1, 0);
}

void pulse_mod_trigger(uint16_t osc) {
    //float mod_sr = (float)AMY_SAMPLE_RATE / (float)AMY_BLOCK_SIZE;
    //float freq = freq_of_logfreq(synth[osc]->logfreq);
    //float period = 1. / (freq/mod_sr);
    //synth[osc]->step = period * synth[osc]->phase;
}

// dpwe sez to use this method for low-freq mod pulse still 
SAMPLE compute_mod_pulse(uint16_t osc) {
    // do BW pulse gen at SR=44100/64
    SAMPLE sample;
    if(msynth[osc]->duty < 0.001f || msynth[osc]->duty > 0.999) msynth[osc]->duty = 0.5;
    if(synth[osc]->phase >= F2P(msynth[osc]->duty)) {
        sample = F2S(1.0f);
    } else {
        sample = F2S(-1.0f);
    }
    float mod_sr = (float)AMY_SAMPLE_RATE / (float)AMY_BLOCK_SIZE;  // samples per sec / samples per call = calls per sec
    float freq = freq_of_logfreq(msynth[osc]->logfreq);
    synth[osc]->phase = P_WRAPPED_SUM(synth[osc]->phase, F2P(freq / mod_sr));  // cycles per sec / calls per sec = cycles per call
    return MULA_SS(sample, F2S(msynth[osc]->amp));
}


/* Saw waves */
void saw_note_on(uint16_t osc, int8_t direction_notused, float freq) {
    // Dummy, now done lazily.
    synth[osc]->lut = NULL;
}

void _saw_note_on(uint16_t osc) {
    //printf("saw_note_on: time %lld osc %d freq %f logfreq %f amp %f last_amp %f phase %f\n", amy_global.total_blocks*AMY_BLOCK_SIZE, osc, freq, synth[osc]->logfreq, msynth[osc]->amp, msynth[osc]->last_amp, P2F(synth[osc]->phase));
    if (synth[osc]->lut == NULL) {
        float freq = freq_of_logfreq(msynth[osc]->logfreq);
        float period_samples = ((float)AMY_SAMPLE_RATE / freq);
        synth[osc]->lut = choose_from_lutset(period_samples, saw_fxpt_lutset);
    }
}

void saw_down_note_on(uint16_t osc, float freq) {
    saw_note_on(osc, -1, freq);
}
void saw_up_note_on(uint16_t osc, float freq) {
    saw_note_on(osc, 1, freq);
}

SAMPLE render_saw(SAMPLE* buf, uint16_t osc, int8_t direction) {
    _saw_note_on(osc);
    return render_lpf_lut(buf, osc, false, direction, /* dc offset */ 0);
    //printf("render_saw: time %lld osc %d buf[]=%f %f %f %f %f %f %f %f\n",
    //       amy_global.total_blocks*AMY_BLOCK_SIZE, osc, S2F(buf[0]), S2F(buf[1]), S2F(buf[2]), S2F(buf[3]), S2F(buf[4]), S2F(buf[5]), S2F(buf[6]), S2F(buf[7]));
}

SAMPLE render_saw_down(SAMPLE* buf, uint16_t osc) {
    return render_saw(buf, osc, -1);
}
SAMPLE render_saw_up(SAMPLE* buf, uint16_t osc) {
    return render_saw(buf, osc, 1);
}


void saw_mod_trigger(uint16_t osc) {
    //float mod_sr = (float)AMY_SAMPLE_RATE / (float)AMY_BLOCK_SIZE;
    //float freq = freq_of_logfreq(synth[osc]->logfreq);
    //float period = 1. / (freq/mod_sr);
    //synth[osc]->step = period * synth[osc]->phase;
}

void saw_up_mod_trigger(uint16_t osc) {
    saw_mod_trigger(osc);
}
void saw_down_mod_trigger(uint16_t osc) {
    saw_mod_trigger(osc);
}

// TODO -- this should use dpwe code
SAMPLE compute_mod_saw(uint16_t osc, int8_t direction) {
    // Saw waveform is just the phasor.
    SAMPLE sample = SHIFTL(P2S(synth[osc]->phase), 1) - F2S(1.0f);
    float mod_sr = (float)AMY_SAMPLE_RATE / (float)AMY_BLOCK_SIZE;  // samples per sec / samples per call = calls per sec
    float freq = freq_of_logfreq(msynth[osc]->logfreq);
    synth[osc]->phase = P_WRAPPED_SUM(synth[osc]->phase, F2P(freq / mod_sr));  // cycles per sec / calls per sec = cycles per call
    return MULA_SS(sample, direction * F2S(msynth[osc]->amp));
}

SAMPLE compute_mod_saw_down(uint16_t osc) {
    return compute_mod_saw(osc, -1);
}

SAMPLE compute_mod_saw_up(uint16_t osc) {
    return compute_mod_saw(osc, 1);
}



/* triangle wave */
void triangle_note_on(uint16_t osc, float freq) {
    // switched to lazy
    synth[osc]->lut = NULL;
}
    
void _triangle_note_on(uint16_t osc, float freq) {
    if (synth[osc]->lut == NULL) {
        float period_samples = (float)AMY_SAMPLE_RATE / freq;
        synth[osc]->lut = choose_from_lutset(period_samples, triangle_fxpt_lutset);
    }
}

SAMPLE render_triangle(SAMPLE* buf, uint16_t osc) {
    float freq = freq_of_logfreq(msynth[osc]->logfreq);
    _triangle_note_on(osc, freq);
    PHASOR step = F2P(freq / (float)AMY_SAMPLE_RATE);  // cycles per sec / samples per sec -> cycles per sample
    SAMPLE amp = F2S(msynth[osc]->amp);
    SAMPLE last_amp = F2S(msynth[osc]->last_amp);
    SAMPLE max_value;
    synth[osc]->phase = render_lut(buf, synth[osc]->phase, step, last_amp, amp, synth[osc]->lut, &max_value);
    msynth[osc]->last_amp = msynth[osc]->amp;
    return max_value;
}

void triangle_mod_trigger(uint16_t osc) {
    // float mod_sr = (float)AMY_SAMPLE_RATE / (float)AMY_BLOCK_SIZE;
    // float freq = freq_of_logfreq(synth[osc]->logfreq);
    // float period = 1. / (freq/mod_sr);
    // synth[osc]->step = period * synth[osc]->phase;
}

// TODO -- this should use dpwe code 
SAMPLE compute_mod_triangle(uint16_t osc) {
    // Offset phase by 1/4 cycle for Triangle waveform in "sine phase" (starts at 0).
    SAMPLE sample = SHIFTL(P2S(synth[osc]->phase), 2) + F2S(1.0f);  // 1..5
    if (sample > F2S(4.0f))  sample -= F2S(4.0f);  // 1..4/0..1
    if (sample > F2S(2.0f))  sample = F2S(4.0f) - sample;  // 0..2..0
    sample -= F2S(1.0f);  // -1 .. 1
    float mod_sr = (float)AMY_SAMPLE_RATE / (float)AMY_BLOCK_SIZE;  // samples per sec / samples per call = calls per sec
    float freq = freq_of_logfreq(msynth[osc]->logfreq);
    synth[osc]->phase = P_WRAPPED_SUM(synth[osc]->phase, F2P(freq / mod_sr));  // cycles per sec / calls per sec = cycles per call
    return MULA_SS(sample, F2S(msynth[osc]->amp));
}


/* FM */
// NB this uses new lingo for step, skip, phase etc
void fm_sine_note_on(uint16_t osc, uint16_t algo_osc) {
    // lazy
    synth[osc]->lut = NULL;
}

void _fm_sine_note_on(uint16_t osc, float freq) {
    if (synth[osc]->lut == NULL) {
        float period_samples = (float)AMY_SAMPLE_RATE / freq;
        synth[osc]->lut = choose_from_lutset(period_samples, sine_fxpt_lutset);
    }
}

SAMPLE render_fm_sine(SAMPLE* buf, uint16_t osc, SAMPLE* mod, SAMPLE feedback_level, uint16_t algo_osc, SAMPLE mod_amp) {
    if(AMY_IS_SET(synth[osc]->logratio)) {
        msynth[osc]->logfreq = msynth[algo_osc]->logfreq + synth[osc]->logratio;
    }
    float freq = freq_of_logfreq(msynth[osc]->logfreq);
    _fm_sine_note_on(osc, freq);
    PHASOR step = F2P(freq / (float)AMY_SAMPLE_RATE);  // cycles per sec / samples per sec -> cycles per sample
    SAMPLE amp = MUL8_SS(F2S(msynth[osc]->amp), mod_amp);
    SAMPLE last_amp = MUL8_SS(F2S(msynth[osc]->last_amp), mod_amp);
    SAMPLE max_value;
    if (feedback_level > 0 && mod)
        synth[osc]->phase = render_lut_fm_fb(buf, synth[osc]->phase, step,
                                            last_amp, amp,
                                            synth[osc]->lut,
                                            mod, feedback_level, synth[osc]->last_two, &max_value);
    else if (feedback_level > 0)
        synth[osc]->phase = render_lut_fb(buf, synth[osc]->phase, step,
                                         last_amp, amp,
                                         synth[osc]->lut,
                                         feedback_level, synth[osc]->last_two, &max_value);
    else if (mod)
        synth[osc]->phase = render_lut_fm(buf, synth[osc]->phase, step,
                                         last_amp, amp,
                                         synth[osc]->lut,
                                         mod, &max_value);
    else
        synth[osc]->phase = render_lut(buf, synth[osc]->phase, step,
                                      last_amp, amp,
                                      synth[osc]->lut, &max_value);

    msynth[osc]->last_amp = msynth[osc]->amp;
    return max_value;
}

/* sine */
void sine_note_on(uint16_t osc, float freq) {
    // functionality now allocated lazily
    synth[osc]->lut = NULL;
}

void _sine_note_on(uint16_t osc, float freq) {
    //fprintf(stderr, "sine_note_on: time %f osc %d freq %f\n", amy_global.total_blocks*AMY_BLOCK_SIZE / (float)AMY_SAMPLE_RATE, osc, freq_of_logfreq(synth[osc]->logfreq_coefs[0]));
    // There's really only one sine table, but for symmetry with the other ones...
    if (synth[osc]->lut == NULL) {
        float period_samples = (float)AMY_SAMPLE_RATE / freq;
        synth[osc]->lut = choose_from_lutset(period_samples, sine_fxpt_lutset);
    }
}

SAMPLE render_sine(SAMPLE* buf, uint16_t osc) { 
    float freq = freq_of_logfreq(msynth[osc]->logfreq);
    _sine_note_on(osc, freq);
    PHASOR step = F2P(freq / (float)AMY_SAMPLE_RATE);  // cycles per sec / samples per sec -> cycles per sample
    SAMPLE amp = F2S(msynth[osc]->amp);
    SAMPLE last_amp = F2S(msynth[osc]->last_amp);
    //fprintf(stderr, "render_sine: time %f osc %d freq %f last_amp %f amp %f\n", amy_global.total_blocks*AMY_BLOCK_SIZE / (float)AMY_SAMPLE_RATE, osc, AMY_SAMPLE_RATE * P2F(step), S2F(last_amp), S2F(amp));
    SAMPLE max_value;
    synth[osc]->phase = render_lut(buf, synth[osc]->phase, step, last_amp, amp, synth[osc]->lut, &max_value);
    msynth[osc]->last_amp = msynth[osc]->amp;
    return max_value;
}


// TOOD -- not needed anymore
SAMPLE compute_mod_sine(uint16_t osc) { 
    // One sample pulled out of render_lut.
    float freq = freq_of_logfreq(msynth[osc]->logfreq);
    _sine_note_on(osc, freq);
    const LUT *lut = synth[osc]->lut;
    if (lut == NULL) return 0;   // Avoid bus error if somehow the osc is not set up yet.
    int lut_mask = lut->table_size - 1;
    int lut_bits = lut->log_2_table_size;
    int16_t base_index = INT_OF_P(synth[osc]->phase, lut_bits);
    SAMPLE frac = S_FRAC_OF_P(synth[osc]->phase, lut_bits);
    LUTSAMPLE b = lut->table[base_index];
    LUTSAMPLE c = lut->table[(base_index + 1) & lut_mask];
    SAMPLE sample = L2S(b) + MUL0_SS(L2S(c - b), frac);
    float mod_sr = (float)AMY_SAMPLE_RATE / (float)AMY_BLOCK_SIZE;  // samples per sec / samples per call = calls per sec
    synth[osc]->phase = P_WRAPPED_SUM(synth[osc]->phase, F2P(freq / mod_sr));  // cycles per sec / calls per sec = cycles per call
    return MULA_SS(sample, F2S(msynth[osc]->amp));
}

void sine_mod_trigger(uint16_t osc) {
    sine_note_on(osc, freq_of_logfreq(msynth[osc]->logfreq));
}


// On the RP2040, rand and mrand48 etc. don't work on core1.
// I think this is to do with the implicit state - it's not repeatable, nor strictly safe,
// to potentially have two cores accessing the same rand_state memory at the same time.
// However, we just want random values, and any bit pattern is acceptable.
// So we make our own implementation of mrand48, and don't worry about it being called from both cores.

static uint64_t rand_state = 0;
static const uint64_t a = 0x5deece66dL;
static const uint32_t c = 0xb;

void my_srand48(uint32_t seedval) {
    rand_state = ((uint64_t)seedval) << 32L;
}

static inline int32_t my_mrand48(void) {
    // per https://www.ibm.com/docs/en/zos/2.4.0?topic=functions-mrand48-pseudo-random-number-generator
    rand_state = (a * rand_state + c) & 0x0000ffffffffffffL;
    return (int32_t)rand_state;
}

#if (defined PICO_RP2350) || (defined PICO_RP2040)
#include "pico/rand.h"
#endif

// Returns a SAMPLE between -1 and 1.
inline static SAMPLE amy_get_random() {
#ifndef AMY_USE_FIXEDPOINT
    return ((float)drand48()) - 0.5f;
#else
    //assert(RAND_MAX == 2147483647); // 2^31 - 1
    return SHIFTR((SAMPLE)my_mrand48(), (32 - S_FRAC_BITS));  // - F2S(0.5f);
#endif
}

/* noise */

void noise_note_on(uint16_t osc) {
    synth[osc]->last_two[0] = 0;
    synth[osc]->last_two[1] = 0;
}

SAMPLE render_noise(SAMPLE *buf, uint16_t osc) {
    SAMPLE amp = F2S(msynth[osc]->amp);
    SAMPLE max_value = 0;
    SAMPLE last_white = synth[osc]->last_two[0];
    SAMPLE last_last_white = synth[osc]->last_two[1];
    for(uint16_t i=0;i<AMY_BLOCK_SIZE;i++) {
        SAMPLE white = MULA_SS(amy_get_random(), amp);
        // Two-zero LPF to make the noise a little more pink, closer to Juno noise.
        SAMPLE value = white + last_white + last_white + last_last_white;
        last_last_white = last_white;
        last_white = white;
        buf[i] += value;
        if (value < 0) value = -value;
        if (value > max_value) max_value = value;
    }
    synth[osc]->last_two[0] = last_white;
    synth[osc]->last_two[1] = last_last_white;
    return max_value;
}

SAMPLE compute_mod_noise(uint16_t osc) {
    float mod_sr = (float)AMY_SAMPLE_RATE / (float)AMY_BLOCK_SIZE;
    float freq = freq_of_logfreq(msynth[osc]->logfreq);
    float fstep = freq / mod_sr;
    SAMPLE amp = F2S(msynth[osc]->amp);
    PHASOR starting_phase = synth[osc]->phase;
    synth[osc]->phase = P_WRAPPED_SUM(synth[osc]->phase, F2P(fstep));  // cycles per sec / calls per sec = cycles per call
    if (fstep > 1.0f || synth[osc]->phase < starting_phase) {
        // phase wrapped, take new sample.
        synth[osc]->last_two[0] = MULA_SS(amy_get_random(), amp);
    }
    //printf("mod_noise: time %lld fstep %f samp %f\n", amy_global.total_blocks*AMY_BLOCK_SIZE, fstep, S2F(synth[osc]->last_two[0]));
    return synth[osc]->last_two[0];
}

/* silent, i.e. just apply envelope */
SAMPLE render_envelope(SAMPLE *buf, uint16_t osc) {
    SAMPLE incoming_amp = F2S(msynth[osc]->last_amp);
    SAMPLE ending_amp = F2S(msynth[osc]->amp);
    SAMPLE max_value = 0;
    SAMPLE current_amp = incoming_amp;
    SAMPLE incremental_amp = SHIFTR(ending_amp - incoming_amp, BLOCK_SIZE_BITS);
    for(uint16_t i = 0; i < AMY_BLOCK_SIZE; i++) {
        SAMPLE value = MULA_SS(buf[i], current_amp);
        buf[i] = value;
        if (value < 0) value = -value;
        if (value > max_value) max_value = value;
        current_amp += incremental_amp;
    }
    msynth[osc]->last_amp = msynth[osc]->amp;
    return max_value;
}


/* partial */

void partial_note_on(uint16_t osc) {
    synth[osc]->lut = NULL;
}

void _partial_note_on(uint16_t osc, float freq) {
    if (synth[osc]->lut == NULL) {
        float period_samples = (float)AMY_SAMPLE_RATE / freq;
        synth[osc]->lut = choose_from_lutset(period_samples, sine_fxpt_lutset);
    }
}

SAMPLE render_partial(SAMPLE * buf, uint16_t osc) {
    float freq = freq_of_logfreq(msynth[osc]->logfreq);
    _partial_note_on(osc, freq);
    PHASOR step = F2P(freq / (float)AMY_SAMPLE_RATE);  // cycles per sec / samples per sec -> cycles per sample
    SAMPLE amp = F2S(msynth[osc]->amp);
    SAMPLE last_amp = F2S(msynth[osc]->last_amp);
    //printf("render_partial: time %.3f logfreq %f freq %f last_amp %f amp %f step %f\n", (float)amy_global.total_blocks*AMY_BLOCK_SIZE/(float)AMY_SAMPLE_RATE, msynth[osc]->logfreq, freq, S2F(last_amp), S2F(amp), P2F(step) * synth[osc]->lut->table_size);
    SAMPLE max_value;
    synth[osc]->phase = render_lut(buf, synth[osc]->phase, step, last_amp, amp, synth[osc]->lut, &max_value);
    msynth[osc]->last_amp = msynth[osc]->amp;
    return max_value;
}

void partial_note_off(uint16_t osc) {
    synth[osc]->substep = 2;
    AMY_UNSET(synth[osc]->note_on_clock);
    synth[osc]->note_off_clock = amy_global.total_blocks*AMY_BLOCK_SIZE;
    msynth[osc]->last_amp = 0;
    synth[osc]->status = SYNTH_OFF;
}


#define MAX_KS_BUFFER_LEN 802 // 44100/55  -- 55Hz (A1) lowest we can go for KS
SAMPLE ** ks_buffer;
uint8_t ks_polyphony_index;


/* karplus-strong */

SAMPLE render_ks(SAMPLE * buf, uint16_t osc) {
    SAMPLE half = MUL0_SS(F2S(0.5f), F2S(synth[osc]->feedback));
    SAMPLE amp = F2S(msynth[osc]->amp);
    float freq = freq_of_logfreq(msynth[osc]->logfreq);
    SAMPLE max_value = 0;
    if(freq >= 55) { // lowest note we can play
        uint16_t buflen = (uint16_t)(AMY_SAMPLE_RATE / freq);
        for(uint16_t i = 0; i < AMY_BLOCK_SIZE; i++) {
            uint16_t index = (uint16_t)(synth[osc]->step);
            SAMPLE sample = ks_buffer[ks_polyphony_index][index];
            ks_buffer[ks_polyphony_index][index] =                 
                SMULR7(
                    (ks_buffer[ks_polyphony_index][index] + ks_buffer[ks_polyphony_index][(index + 1) % buflen]),
                    half);
            synth[osc]->step = (index + 1) % buflen;
            SAMPLE value = SMULR7(sample, amp);
            buf[i] += value;
            if (i == 0) {
                max_value = value;
            } else {
                if (value > max_value) max_value = value;
                else if (-value > max_value) max_value = -value;
            }
        }
    }
    //fprintf(stderr, "render_ks time %u osc %d freq %.1f amp %.3f maxval %.3f\n", amy_global.total_blocks*AMY_BLOCK_SIZE, osc, freq, S2F(amp), S2F(max_value));
    return max_value;
}

void ks_note_on(uint16_t osc) {
    float freq = freq_of_logfreq(msynth[osc]->logfreq);
    if(freq <= 1.f) freq = 1.f;
    uint16_t buflen = (uint16_t)(AMY_SAMPLE_RATE / freq);
    if(buflen > MAX_KS_BUFFER_LEN) buflen = MAX_KS_BUFFER_LEN;
    // init KS buffer with noise up to max
    SAMPLE sum = 0;
    for(uint16_t i = 0; i < buflen; i++) {
        SAMPLE val = amy_get_random();
        ks_buffer[ks_polyphony_index][i] = val;
        sum += val;
    }
    // Remove dc, to avoid ending up with a dc-offset residual.
    SAMPLE mean = sum / buflen;
    for(uint16_t i = 0; i < buflen; i++) {
        ks_buffer[ks_polyphony_index][i] -= mean;
    }
    ks_polyphony_index++;
    if(ks_polyphony_index == AMY_KS_OSCS) ks_polyphony_index = 0;
    //fprintf(stderr, "ks_note_on: osc %d buflen %d poly_index %d\n", osc, buflen, ks_polyphony_index);
}

void ks_note_off(uint16_t osc) {
    msynth[osc]->amp = 0;
}


void ks_init(void) {
    // 6ms buffer
    ks_polyphony_index = 0;
    ks_buffer = (SAMPLE**) malloc(sizeof(SAMPLE*)*AMY_KS_OSCS);
    for(int i=0;i<AMY_KS_OSCS;i++) ks_buffer[i] = (SAMPLE*)malloc(sizeof(float)*MAX_KS_BUFFER_LEN); 
}

void ks_deinit(void) {
    for(int i=0;i<AMY_KS_OSCS;i++) free(ks_buffer[i]);
    free(ks_buffer);
}

// --------- wavetable ----------

#ifdef AMY_WAVETABLE
void wavetable_note_on(uint16_t osc, float freq) {
    //fprintf(stderr, "wavetable_note_on: time %f osc %d freq %f\n", amy_global.total_blocks*AMY_BLOCK_SIZE / (float)AMY_SAMPLE_RATE, osc, freq_of_logfreq(synth[osc]->logfreq_coefs[0]));
    //float period_samples = (float)AMY_SAMPLE_RATE / freq;
    //synth[osc]->lut = wavetable_lut;  // TODO(dpwe): choose based on synth[osc]->preset.
}

// Structure of waveeditonlie wavetables.
const int CYCLES_PER_WAVETABLE = 64;
const int WAVETABLE_SAMPLES_PER_CYCLE = 256;
const int WAVETABLE_LOG2_SAMPLES_PER_CYCLE = 8;

SAMPLE render_wavetable(SAMPLE* buf, uint16_t osc) {
    float freq = freq_of_logfreq(msynth[osc]->logfreq);
    PHASOR step = F2P(freq / (float)AMY_SAMPLE_RATE);  // cycles per sec / samples per sec -> cycles per sample
    SAMPLE amp = F2S(msynth[osc]->amp);
    SAMPLE last_amp = F2S(msynth[osc]->last_amp);
    //fprintf(stderr, "render_wavetable: time %f osc %d freq %f last_amp %f amp %f preset %d\n", amy_global.total_blocks*AMY_BLOCK_SIZE / (float)AMY_SAMPLE_RATE, osc, AMY_SAMPLE_RATE * P2F(step), S2F(last_amp), S2F(amp), synth[osc]->preset);
    SAMPLE max_value;
    float interp = MAX(0, MIN(CYCLES_PER_WAVETABLE - 1, (CYCLES_PER_WAVETABLE - 1) * msynth[osc]->duty));  // Don't try to interp beyond end of table.  An N-waveform table can be interpolated from 0 to (N-1-eps).
    int table = MIN((int)floor(interp), CYCLES_PER_WAVETABLE - 2);  // always need both this wavetable and the next one.
    interp = interp - table;  // fractional part, normally < 1.0, but == 1.0 for very end of table.
    int wavetable_preset = AMY_IS_SET(synth[osc]->preset)
        ? (int)synth[osc]->preset
        : (int)pcm_wavetable_base;
    int wavetable_samples_per_table = (pcm_wavetable_len > 0) ? (int)pcm_wavetable_len : 16384;
    uint32_t sample_length = 0;
    const int16_t *wavetable_sample_ram =
        pcm_get_sample_ram_for_preset((uint16_t)wavetable_preset, &sample_length);
    if ((wavetable_sample_ram == NULL || sample_length < (uint32_t)wavetable_samples_per_table) &&
        wavetable_preset != (int)pcm_wavetable_base) {
        wavetable_sample_ram = pcm_get_sample_ram_for_preset(pcm_wavetable_base, &sample_length);
    }
    if (wavetable_sample_ram == NULL || sample_length < (uint32_t)wavetable_samples_per_table) {
        return 0;
    }
    LUT wavetable_lut = {wavetable_sample_ram, WAVETABLE_SAMPLES_PER_CYCLE, WAVETABLE_LOG2_SAMPLES_PER_CYCLE, 0, 1.0f};
    wavetable_lut.table += table * WAVETABLE_SAMPLES_PER_CYCLE;
    // don't update phase in the first call to render_lut, so second call uses the same phase
    SAMPLE interp_a = F2S(1.0f - interp);
    SAMPLE interp_b = F2S(interp);
    // If we used last_duty, we could actually smoothly interpolate the waveshape crossfade too (except across table boundaries).
    render_lut(buf, synth[osc]->phase, step, SMULR7(last_amp, interp_a), SMULR7(amp, interp_a), &wavetable_lut, &max_value);
    // Point to next cycle.
    wavetable_lut.table += WAVETABLE_SAMPLES_PER_CYCLE;
    synth[osc]->phase = render_lut(buf, synth[osc]->phase, step, SMULR7(last_amp, interp_b), SMULR7(amp, interp_b), &wavetable_lut, &max_value);
    msynth[osc]->last_amp = msynth[osc]->amp;
    return max_value;
}
#endif
