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 [email protected].
To post to this group, send email to [email protected].
Visit this group at https://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.
