mpy.polys.numberfields.galoisgroups import Resolvent
>>> X = symbols('X0 X1 X2 X3')
>>> F = X[0]*X[2] + X[1]*X[3]
>>> s = [Permutation([0, 1, 2, 3]), Permutation([1, 0, 2, 3]),
... Permutation([3, 1, 2, 0])]
>>> R = Resolvent(F, X, s)

This resolvent has three roots, which are the conjugates of ``F`` under the
three permutations in ``s``:

>>> R.root_lambdas[0](*X)
X0*X2 + X1*X3
>>> R.root_lambdas[1](*X)
X0*X3 + X1*X2
>>> R.root_lambdas[2](*X)
X0*X1 + X2*X3

Resolvents are useful for computing Galois groups. Given a polynomial $T$
of degree $n$, we will use a resolvent $R$ where $Gal(T) \leq G \leq S_n$.
We will then want to substitute the roots of $T$ for the variables $X_i$
in $R$, and study things like the discriminant of $R$, and the way $R$
factors over $\mathbb{Q}$.

From the symmetry in $R$'s construction, and since $Gal(T) \leq G$, we know
from Galois theory that the coefficients of $R$ must lie in $\mathbb{Z}$.
This allows us to compute the coefficients of $R$ by approximating the
roots of $T$ to sufficient precision, plugging these values in for the
variables $X_i$ in the coefficient expressions of $R$, and then simply
rounding to the nearest integer.

In order to determine a sufficient precision for the roots of $T$, this
``Resolvent`` class imposes certain requirements on the form ``F``. It
could be possible to design a different ``Resolvent`` class, that made
different precision estimates, and different assumptions about ``F``.

``F`` must be homogeneous, and all terms must have unit coefficient.
Furthermore, if $r$ is the number of terms in ``F``, and $t$ the total
degree, and if $m$ is the number of conjugates of ``F``, i.e. the number
of permutations in ``s``, then we require that $m < r 2^t$. Again, it is
not impossible to work with forms ``F`` that violate these assumptions, but
this ``Resolvent`` class requires them.

Since determining the integer coefficients of the resolvent for a given
polynomial $T$ is one of the main problems this class solves, we take some
time to explain the precision bounds it uses.

The general problem is:
Given a multivariate polynomial $P \in \mathbb{Z}[X_1, \ldots, X_n]$, and a
bound $M \in \mathbb{R}_+$, compute an $\varepsilon > 0$ such that for any
complex numbers $a_1, \ldots, a_n$ with $|a_i| < M$, if the $a_i$ are
approximated to within an accuracy of $\varepsilon$ by $b_i$, that is,
$|a_i - b_i| < \varepsilon$ for $i = 1, \ldots, n$, then
$|P(a_1, \ldots, a_n) - P(b_1, \ldots, b_n)| < 1/2$. In other words, if it
is known that $P(a_1, \ldots, a_n) = c$ for some $c \in \mathbb{Z}$, then
$P(b_1, \ldots, b_n)$ can be rounded to the nearest integer in order to
determine $c$.

To derive our error bound, consider the monomial $xyz$. Defining
$d_i = b_i - a_i$, our error is
$|(a_1 + d_1)(a_2 + d_2)(a_3 + d_3) - a_1 a_2 a_3|$, which is bounded
above by $|(M + \varepsilon)^3 - M^3|$. Passing to a general monomial of
total degree $t$, this expression is bounded by
$M^{t-1}\varepsilon(t + 2^t\varepsilon/M)$ provided $\varepsilon < M$,
and by $(t+1)M^{t-1}\varepsilon$ provided $\varepsilon < M/2^t$.
But since our goal is to make the error less than $1/2$, we will choose
$\varepsilon < 1/(2(t+1)M^{t-1})$, which implies the condition that
$\varepsilon < M/2^t$, as long as $M \geq 2$.

Passing from the general monomial to the general polynomial is easy, by
scaling and summing error bounds.

In our specific case, we are given a homogeneous polynomial $F$ of
$r$ terms and total degree $t$, all of whose coefficients are $\pm 1$. We
are given the $m$ permutations that make the conjugates of $F$, and
we want to bound the error in the coefficients of the monic polynomial
$R(Y)$ having $F$ and its conjugates as roots (i.e. the resolvent).

For $j$ from $1$ to $m$, the coefficient of $Y^{m-j}$ in $R(Y)$ is the
$j$th elementary symmetric polynomial in the conjugates of $F$. This sums
the products of these conjugates, taken $j$ at a time, in all possible
combinations. There are $\binom{m}{j}$ such combinations, and each product
of $j$ conjugates of $F$ expands to a sum of $r^j$ terms, each of unit
coefficient, and total degree $jt$. An error bound for the $j$th coeff of
$R$ is therefore
$$\binom{m}{j} r^j (jt + 1) M^{jt - 1} \varepsilon$$
When our goal is to evaluate all the coefficients of $R$, we will want to
use the maximum of these error bounds. It is clear that this bound is
strictly increasing for $j$ up to the ceiling of $m/2$. After that point,
the first factor $\binom{m}{j}$ begins to decrease, while the others
continue to increase. However, the binomial coefficient never falls by more
than a factor of $1/m$ at a time, so our assumptions that $M \geq 2$ and
$m < r 2^t$ are enough to tell us that the constant coefficient of $R$,
i.e. that where $j = m$, has the largest error bound. Therefore we can use
$$r^m (mt + 1) M^{mt - 1} \varepsilon$$
as our error bound for all the coefficients.

Note that this bound is also (more than) adequate to determine whether any
of the roots of $R$ is an integer. Each of these roots is a single
conjugate of $F$, which contains less error than the trace, i.e. the
coefficient of $Y^{m - 1}$. By rounding the roots of $R$ to the nearest
integers, we therefore get all the candidates for integer roots of $R$. By
plugging these candidates into $R$, we can check whether any of them
actually is a root.

Note: We take the definition of resolvent from Cohen, but the error bound
is ours.

References
==========

.. [1] Cohen, H. *A Course in Computational Algebraic Number Theory*.
   (Def 6.3.2)

c