"""
This is an implementation of the algorithm described in
Lehobey's
"Resolvent Computations by Resultants Without Extraneous Powers"
"""


def interpolating_functions(f):
  """
  Build a list of interpolating functions for f

  We require f : K[x1,...,xn][x]
  """
  R = f.parent().base_ring()
  xs = R.gens()

  n = f.degree()

  out = [f]
  for k in range(f.degree()):
    fk = (out[-1] - out[-1].subs(x=xs[n-k-1]))/(x - xs[n-k-1])
    out += [f.parent()(fk.simplify_full())]

  return out[::-1]


def stabilizer(G,p):
  """
  Compute the stabilizer of p by G
  """
  elems = []
  for g in G:
    if p * g == p:
      elems += [g]
  return G.subgroup(elems)



def truncated_root(p,r,d):
  """
  Compute q so that q^r = p (working mod x^d)

  Assumes p : A[t] has constant term 1
  and that such a q : A[t] actually exists!
  """

  r = int(r)

  t = p.variables()[0]

  n = p.degree()
  # print("n: ", n)

  p = p.truncate(d)

  # print("truncating root -- ")
  # print("p: ", p)
  # print("r: ", r)
  # print("d: ", d)

  ps = p.coefficients(sparse=False)

  # these will be the coefficients of q
  qs = [1]

  for k in range(n // r):
    qs += [1/(k+1) * sum([(k+1 - (r+1)*j) * qs[j] * ps[k+1-j] / r for j in range(k+1)])]

  return sum([qs[j] * t^j for j in range(len(qs))])



def resolvent(f,Theta):
  """
  Compute mathcal{L}_{Theta,f} as per the paper

  Assumes f : K[x1,...,xn][x] and Theta : K[x1,...,xn]
  """
  R = Theta.parent()
  K = R.base_ring()
  xs = R.gens()

  SIterated.<t> = PolynomialRing(R)
  S = SIterated.flattening_morphism().codomain()

  T = K[t]

  Rj = (t - SIterated(Theta)).reverse()
  Rj = S(Rj)

  fs = interpolating_functions(f)

  HprevOrder = 1

  n = f.degree()

  for j in range(1,n+1):
    print("j: ", j)
    Sj = SymmetricGroup(j)
    Hj = stabilizer(Sj,Theta)
    dj = factorial(j) / Hj.order()
    mj = Hj.order() / HprevOrder

    # print("Sj: ", Sj)
    # print("Hj: ", Hj)
    # print("dj: ", dj)
    # print("mj: ", mj)

    # update the previous order for the next cycle
    HprevOrder = Hj.order()

    # there's an annoying off-by-one error with the variable names
    # compared to everything else
    fj = S(fs[j].subs(x=xs[j-1]))

    # print("fj: ", fj)

    res = fj.resultant(Rj, S(xs[j-1]))

    # print("res: ", res)

    Rj = truncated_root(SIterated(res),mj,dj+1)
    # print("Rj: ", Rj)

    Rj = S(Rj)

    # print("-----")

  # print("about to return: ")
  # print(Rj)
  return T(Rj).reverse()


def solveByRadicals(f):
  """
  Compute a root of f using radicals
  """

  n = int(f.degree(x))
  K.<w> = CyclotomicField(n)

  R = PolynomialRing(K,n,'x')
  xs = R.gens()

  R1 = R[x]
  f = R1(f)

  Theta = sum(xs[k] * w^(k) for k in range(n))
  print("Theta: ", Theta)

  # Theta^5 is preserved under the action of the galois group,
  # while Theta itself is an eigenvector with eigenvalue w
  L = resolvent(f,Theta^n)

  print("L: ", L)

  psis = [-list(f)[-2]] + L.roots(multiplicities=False)

  print("psis: ")
  print([psi.n() for psi in psis])

  thetas = [psi^(1/n) for psi in psis]

  r = sum(thetas) / n

  print("purported root r: ", r)
  print("f(r): ", f(r).n())


R = QQ[x]


# f = x^3 + x^2 - 2*x - 1
# f = x^5 + x^4 -4*x^3 - 3*x^2 + 3*x + 1

# this works for degree 3 polynomials
# fs = [x^3 - x^2 - a*x + b for (a,b) in [(26,-41), (32,79), (34,61), (36,4), (42,-80), (46,-103)]]

# what about degree 5?
fs = [x^5 + x^4 - 4*x^3 - 3*x^2 + 3*x + 1, # this one doesn't
      x^5 + x^4 - 12*x^3 - 21*x^2 + 1*x + 5, # this one works
      x^5 + x^4 - 16*x^3 + 5*x^2 + 21*x - 9, # this one doesn't work
      x^5 + x^4 - 24*x^3 - 17*x^2 + 41*x - 13]


for f in fs:
  print("f: ", f)
  print("Gf: ", R(f).galois_group().structure_description())
  solveByRadicals(f)
  print("--------------")