// What follows is a quick and dirty implementation of polynomials
// represented as lists of coefficients, and an attempt to implement
// the subresultant pseudo-remainder sequence ("Subresultant PRS")
// according to a combination of Zippel's "Effective Polynomial
// Computation" section 8.6. page 135, and the wikipedia page
// https://en.wikipedia.org/wiki/Polynomial_greatest_common_divisor

function numcmp(a: number, b: number): -1 | 0 | 1 {
  return a < b ? -1 : a === b ? 0 : 1;
}

class poly {
  constructor(public t: number[]) {
    this.t = poly.normalize(this.t);
  }
  toString() {
    if (this.t.length === 0) return '0';
    let rv: string = this.t[0] + '';
    for (let i = 1; i < this.t.length; i++) {
      if (this.t[i] === 0) continue;
      const coe = this.t[i] === 1 ? '' : this.t[i];
      if (i == 1)
        rv = `${coe}x + ${rv}`;
      else
        rv = `${coe}x^${i} + ${rv}`;
    }
    return rv.replace(/\+ -/g, '- ').replace(/ \+ 0$/, '');
  }
  inspect() {
    return this.toString();
  }

  isZero(): boolean {
    // relies on normalization
    return this.t.length === 0;
  }

  // degree of polynomial
  deg(): number {
    // relies on normalization
    return this.t.length - 1;
  }

  // leading coefficient
  lc(): number {
    // relies on normalization
    if (this.t.length === 0) return 0;
    return this.t[this.t.length - 1];
  }

  scale(n: number): poly {
    return new poly(this.t.map(x => x * n));
  }

  static normalize(orig: number[]): number[] {
    const t = [...orig];
    t.forEach((n, i) => {
      if (Math.abs(n) < 1e-10)
        t[i] = 0;
    });
    let zeroes = [...t].reverse().findIndex(x => x !== 0);
    if (zeroes === -1) {
      zeroes = t.length;
    }
    return t.slice(0, t.length - zeroes);
  }

  lead1(): poly {
    return poly.mult(this, new poly([1 / this.lc()]));
  }

  static zero(): poly {
    return new poly([]);
  }

  static mono(mspec: { coeff: number, deg: number }): poly {
    const { coeff, deg } = mspec;
    const rv: number[] = [];
    for (let i = 0; i <= deg; i++) {
      rv[i] = i == deg ? coeff : 0;
    }
    return new poly(rv);
  }

  static add(a: poly, b: poly): poly {
    const rv: number[] = [];
    for (let i = 0; i < a.t.length; i++) {
      rv[i] = (rv[i] || 0) + a.t[i];
    }
    for (let i = 0; i < b.t.length; i++) {
      rv[i] = (rv[i] || 0) + b.t[i];
    }
    return new poly(rv);
  }

  static sub(a: poly, b: poly): poly {
    const rv: number[] = [];
    for (let i = 0; i < a.t.length; i++) {
      rv[i] = (rv[i] || 0) + a.t[i];
    }
    for (let i = 0; i < b.t.length; i++) {
      rv[i] = (rv[i] || 0) - b.t[i];
    }
    return new poly(rv);
  }

  static sgn(a: poly): -1 | 0 | 1 {
    return Math.sign(a.lc()) as (-1 | 0 | 1);
  }

  static cmp(a: poly, b: poly): -1 | 0 | 1 {
    return numcmp(a.deg(), b.deg()) || poly.sgn(poly.sub(a, b));
  }

  private static mult2(x: poly, y: poly): poly {
    const rv: number[] = [];
    x.t.forEach((nx, ix) => {
      y.t.forEach((ny, iy) => {
        const i = ix + iy;
        if (!rv[i]) rv[i] = 0;
        rv[i] += nx * ny;
      });
    });
    return new poly(rv);
  }

  static mult(...ps: poly[]): poly {
    if (ps.length === 0) return pol(1);
    if (ps.length === 1) return ps[0];
    if (ps.length === 2) return poly.mult2(ps[0], ps[1]);
    return poly.mult2(ps[0], poly.mult(...ps.slice(1)));
  }
}

function pol(...t: number[]): poly {
  return new poly(t);
}

function root(t: number) { return pol(-t, 1, 0) }

function div(big: poly, small: poly): { quot: poly, rem: poly } {
  let rem = big;
  let quot = poly.zero();
  let oldld = rem.deg();
  while (rem.deg() >= small.deg()) {
    const mono = poly.mono({
      coeff: rem.lc() / small.lc(),
      deg: rem.deg() - small.deg()
    });
    rem = poly.sub(rem, poly.mult(mono, small));
    quot = poly.add(quot, mono);
    if (!(rem.deg() < oldld))
      throw `leading degree did not decrease during (${big})/(${small}), new rem is ${rem} something went wrong`;
  }
  return { quot, rem };
}

// Compute the subresultant pseudo-remainder sequence.

function subres(p: poly, q: poly): poly[] {
  // r[i] here is F_{i+1} in Zippel
  // d[i] here is δ_i in Zippel
  // gamma[i] here is f_{i+1} in Zippel
  // psi[i] here is h_{i} in Zippel
  // beta[i] here is β_{i+2} in Zippel
  let r: poly[] = [p, q];
  let d: number[] = [];
  let beta: number[] = [];
  let gamma: number[] = [];
  let psi: number[] = [];
  for (let i = 1; r[i].deg() > 0; i++) {
    d[i] = r[i - 1].deg() - r[i].deg();
    gamma[i] = r[i].lc();
    if (i == 1) {
      psi[1] = 1;
      beta[1] = Math.pow(-1, d[1] + 1);
    }
    else {
      // Zippel says: h_i = Math.pow(f_i, δ_{i-1}) / Math.pow(h_{i-1}, d_{i-1} - 1);
      psi[i] = Math.pow(gamma[i - 1], d[i - 1]) / Math.pow(psi[i - 1], d[i - 1] - 1);
      // Zippel says: β_{i} = Math.pow(-1, δ_{i-2} + 1) * f_{i-2} * Math.pow(h_{i-2}, δ_{i-2});
      beta[i] = Math.pow(-1, d[i] + 1) * gamma[i - 1] * Math.pow(psi[i], d[i]);
    }
    r[i + 1] = div(
      r[i - 1].scale(Math.pow(gamma[i], d[i] + 1)),
      r[i]
    ).rem.scale(1 / beta[i]);
  }
  return r;
}

const p = pol(-5, 2, 8, -3, -3, 0, 1, 0, 1);
const q = pol(21, -9, -4, 0, 5, 0, 3);
console.log(subres(p, q));
// I get:
// [ x^8 + x^6 - 3x^4 - 3x^3 + 8x^2 + 2x - 5,
//   3x^6 + 5x^4 - 4x^2 - 9x + 21,
//   15x^4 - 3x^2 + 9,
//   65x^2 + 125.00000000000001x - 245.00000000000003,
//   9326.000000000004x - 12300.000000000005,
//   260708.00000000035 ]