For the problem at hand here, I would work with polynomials in `x` with coefficients in the ring of polynomials in `N` and `M` over the field of algebraic numbers, and do something like the following.
Define the ring of polynomials over the field of algebraic numbers. sage: R.<N, M> = QQbar[] Define the ring of polynomials over the above polynomial ring. sage: S.<x> = R[] Define `a`, `i`, `b`, `c`, `X`, `f` from the question. sage: a = QQbar(3).sqrt() / 2 sage: i = QQbar(I) sage: b = (x - ~x) / (2 * i) sage: c = (a * (x + ~x) + b) / 2 sage: X = (M / 3) * (a + b + c) sage: f = 24 * (X^2 - N * b * c) * x^2 Note that the definition of `b` involves the inverse of `x` (which can be denoted by `~x` or `x^-1`) and therefore lives in the fraction field of `S`, rather than in `S` itself. Therefore, so do `c`, `X`, and `f`. Check if `f` however represents an element in `S` as follows: sage: f in S True and then construct the corresponding element of `S` (we could call it `f` but here we call it `ff`). sage: ff = S(f) Check that `ff` vanishes at `-1`, or equivalently is divisible by `x + 1`. sage: ff(-1) 0 sage: d = x + 1 sage: d.divides(ff) True Perform exact division (to stay in the polynomial ring `S` rather than move to its fraction field). sage: g = ff // d Define `t` as an algebraic number. sage: t = QQbar(exp(-pi*I/3)) Check that `g` vanishes at `t` or equivalently is divisible by `x - t`. sage: g(t) 0 sage: dd = x - t sage: dd.divides(g) True Perform exact division. sage: h = g // dd Inspect the result: sage: h ((-1.000000000000000? - 1.732050807568878?*I)*M^2 + (3.000000000000000? + 5.196152422706632?*I)*N)*x^2 + ((1.000000000000000? - 1.732050807568878?*I)*M^2 + (3.000000000000000? - 5.196152422706632?*I)*N)*x + 2*M^2 + (-6)*N Or as a list (the list of coefficients of `x^j`, for `j` from `0` to the degree): sage: h.list() [2*M^2 + (-6)*N, (1.000000000000000? - 1.732050807568878?*I)*M^2 + (3.000000000000000? - 5.196152422706632?*I)*N, (-1.000000000000000? - 1.732050807568878?*I)*M^2 + (3.000000000000000? + 5.196152422706632?*I)*N] Define somme pretty_print functions to inspect these algebraic numbers: def pretty_print_qqbar(z): a, b = z.real().radical_expression(), z.imag().radical_expression() return "{} + {}*i".format(a, b) def pretty_print_x_coeff(xc): return " + ".join("({})*{}".format(pretty_print_qqbar(xc[m]), m) for m in xc.monomials()) and print `h` in a nice form: print("h = " + "\n + ".join("({})*x^{}".format(pretty_print_x_coeff(xc), k) for k, xc in enumerate(g.list()))) The result is: h = ((-1 + sqrt(3)*i)*M^2 + (3 + -3*sqrt(3)*i)*N)*x^0 + ((3 + sqrt(3)*i)*M^2 + (-3 + 3*sqrt(3)*i)*N)*x^1 + ((3 + -sqrt(3)*i)*M^2 + (-3 + -3*sqrt(3)*i)*N)*x^2 + ((-1 + -sqrt(3)*i)*M^2 + (3 + 3*sqrt(3)*i)*N)*x^3 -- You received this message because you are subscribed to the Google Groups "sage-support" group. To unsubscribe from this group and stop receiving emails from it, send an email to sage-support+unsubscr...@googlegroups.com. To post to this group, send email to sage-support@googlegroups.com. Visit this group at https://groups.google.com/group/sage-support. For more options, visit https://groups.google.com/d/optout.