typedef struct {
  float x_old = 0;
  float Y_old = 0;
  float a = 0;
} waveshaper_state;

/*
 * The function i'm trying to map is: f(x) = sgn(x)*(ax^2)*(ax^2+3)/(ax^2+1)^2
 * Its integral is equal to F(x) = |x|*(ax^2)/(ax^2+1)
 * The function peaks at x = (+/-) sqrt(3/a), and the peak value is (+/-) 9/8
 * I implemented the function in floating point, it's cheaper than fixed point
 *
 * You easily can hack the code. If you do it, please give credit.
 */

float waveshaper_update(float x, waveshaper_state *wss) {
  float x_abs;
  float x_sgn;

  if (x > 0) {
    x_abs = x;
    x_sgn = 1;
  } else {
    x_abs = -x;
    x_sgn = -1;
  }

  float ax2 = (wss->a) * x * x;    // ax^2
  float g_x = 1.0f / (1.0f + ax2); // 1/(b+ax^2)
  float y = x_sgn * ax2 * g_x * g_x * (ax2 + 3.0f);
  float Y = x_abs * ax2 * g_x;

  float dY = Y - (wss->Y_old);
  float dx = x - (wss->x_old);
  float dx_abs = dx > 0 ? dx : -dx;

  float mean = dY / dx;
  (wss->x_old) = x;
  (wss->Y_old) = Y;
  return dx_abs > 0.0001f ? mean : y;
}

void waveshaper_coefs(int32_t squeeze, waveshaper_state *wss) {
  float xmax = 1.05f - (squeeze / 134217728.0f);
  (wss->a) = 3.0f / (xmax * xmax);
}

waveshaper_state satu;

waveshaper_coefs(param_amt, &satu);

outlet_out =
    119304647.0f * waveshaper_update(((float)inlet_in) / 134217728.0f, &satu);
// that 119304647 corresponds (in Q27) to 8/9
// This is because this function peaks at 9/8 before asymptotically going to 1,
// so in order to keep the output in normal range