uint32_t p; // phase
int32_t dp; // delta phase aka freq
int32_t x0, x1, I1_0, I1_1;
int32_t kMI;      // Modulation Index
int32_t aMI, dMI; // interp Modulation Index

// _____________________________________________________________________

// _____________________________________________________________________
p = 0;
dp = 1;
x0 = 0;
x1 = 1;
I1_0 = I1_1 = 0;
aMI = dMI = kMI = 0;
// _____________________________________________________________________

if (!inlet_disable) {
  // _____________________________________________________________________
  MTOFEXTENDED(param_pitch + inlet_pitch, dp);
  dp += ___SMMUL(0x01000000, inlet_freq_space_offset);
  dp >>= 32 - 20 - attr_piwt.LENGTHPOW;
  // (need to calculate x0 - x1 even with over modulation)

  // _____________________________________________________________________
  //                                                     Audio rate loops
  const int rightShiftPM = 9 - attr_piwt.LENGTHPOW;
  const int shiftDisto = attr_piwt.LENGTHPOW - 4;
  for (int j = 0; j < BUFSIZE; j++) {
    if (inlet_disto) {
      p = (inlet_disto_space_offset >> 4) << shiftDisto;
    } else {
      p += dp; // phase increment
      if (p >= 0x40000000) {
        p -= attr_piwt.LENGTH << 20;  // length of one cycle
        x0 -= attr_piwt.LENGTH << 20; // p, xO and x1 must be consistent
      }
    }
    x1 = x0;
    I1_1 = I1_0;

    x0 = p + ((inlet_pm[j] << 3) >> rightShiftPM);

    // differentiation part 1/3
    // we launch the float div here (14 cycles !!)
    // so that integer ops execute in parallel :D
    float inv = ((float)(1 << 27)) / ((float)(x0 - x1));

    // Now we do as many int ops as we can to fill ~14 cycles

    // the interpolation coefficient
    int32_t a = (x0 & 0x000FFFFF) << 11; // 0 -> 1 q31

    // the segment number
    uint32_t i =
        (x0 & (attr_piwt.LENGTHMASK << 20)) >> 20; // [0 attr_piwt.LENGTH-1]
    // polynomial calc of the integrand on segment i at position a
    // the polynomial is
    // I1_0 = Iv + a * v + (a^2/2) dv
    // it is the integrand of v + a*dv

    // a^2/2  q31
    int32_t a2_2 = ___SMMUL(a, a);
    I1_0 = ___SMMUL(attr_piwt.v[i], a);
    I1_0 = ___SMMLA(attr_piwt.v[i + 1] - attr_piwt.v[i], a2_2, I1_0);
    I1_0 = attr_piwt.Iv[i] + (I1_0 << 1);
    int32_t v =
        ___SMMLA(a, attr_piwt.v[i + 1] - attr_piwt.v[i], attr_piwt.v[i] >> 1)
        << (attr_piwt.GAIN + 1);

    // we differentiate the integrand in the interval
    // x0 x1: ie in the interval between now and the preceding
    // input.

    // differentiation part 2/3
    int32_t diff = I1_0 - I1_1;

    // now, the result of our float division should be available.
    // we haven't wasted our time waiting for it.

    // differentiation part 3/3
    if (abs(x0 - x1) > 256) // we avoid ill cases
      outlet_wave[j] = (int32_t)(diff * inv);
    else {
      outlet_wave[j] = v;
    }
  }
} else { // disabled
  if (inlet_disto) {
    for (int j = 0; j < BUFSIZE; j++) {
      outlet_wave[j] = inlet_pm[j];
    }
  } else {
    for (int j = 0; j < BUFSIZE; j++) {
      outlet_wave[j] = 0;
    }
  }
}