Preamble:
I want to calculate the minimal polynomial of a number of the form
$$x=\sum_{i=1}^k \pm a_i^{1/k_i}$$
Where the $a_i$ are algebraic numbers also of this form with a finite expression and $k_i$ are positive integers. I can think of three approaches, all with serious problems:
- Using the resultant repeatedly:
Apply the algorithm recursively until you arrive at a root of a rational number, therefore we can assume that we have the minimal polynomial of all $a_i$. Then compute the minimal polynomial of $\pm a_1^{1/k_1}\pm a_2^{1/k_2}$ using the resultant of their minimal polynomials then repeat for that with $a_3, a_4...a_k$. Let $d_i$ be the degree of the minimal polynomial of $a_i^{1/k_i}$. In general this method requires computing the determinant of square matrices of orders $d_1+d_2$, $d_1+d_2+d_3,\ldots$ up to $d_1+d_2+\ldots+d_k$. According to wikipedia calculating the determinant of a matrix of order $n$ has complexity of the order $O(n^{2.373})$ but practically it's faster to use the LU decomposition method (this is what, for example, the blaze c++ library does) which has a complexity of $O(n^3)$. If we have $d_1=d_2=...=d$ then the complexity is $O((2d)^3+(3d)^3...(kd)^3)=O(d^3k^4)$.
- Using an Integer relation algorithm:
If we know that the degree of $x$ is $n$ then we can run for example the LLL algorithm on ${1,x,x^2...x^n}$. Problems with that include:
Determination of $n$: in the previous example we can get an upper bound of $d^k$ (yikes) so running the LLL algorithm on this would have complexity $O(d^k)$ assuming I'm interpreting wikipedia correctly. If $n$ does not match an upper bound then we either need to run the algorithm again or factor the polynomial.
Precision: if we calculate these values outwards starting at the rational constants (for example $\sqrt[3]{2-\sqrt{2}}\approx\sqrt[3]{0.5858}\approx 0.8367$) we run into the problem that each operation has a tendency to exacerbate the errors of the previous ones. This is particularly problematic when subtracting two large numbers whose difference is small. So we'd need bounds on the error $x^i\pm\epsilon$ where the algorithm still returns the correct result. Edit: since writing this, I found this blog post by Jack Coughlin which details how to estimate this error. This greatly mitigates the problem.
Afterwards we'd need to check if the polynomial found really has $x$ as a root and whether it's irreducible (this is only necessary if we don't know $n$ in advance, then we'd have to use an upper bound).
- Trying some symbolic manipulation, especially as a middle step of another algorithm:
In some cases isolating the term with the highest root might be helpful. For example
$$x=\sqrt{a}+\sqrt[3]{b}$$
$$x^3-3\sqrt a x^2+3ax-a\sqrt a = b$$
$$(x^3+3ax-b)^2=9ax^4-6a^2x^2+a^3$$
I get the impression that this isn't efficient, so maybe we can use a lookup table for small exponents with little nesting.
Applications
Of course I want to implement this in a computer so it does all the number crunching for me. For example, recently I was trying to do a search over graphs composed of triangular lattices so what I did is set the vertices of the graph as points in $\mathbb C$ so I had to check whether several numbers of the form $\eta(a+b\omega)+\mu$ where $w=\frac{1-\sqrt{-3}}2$ are equal to each other with sympy but since it stores expressions "as is" there's no trivial way to check equality, so I needed to check if the minimal polynomial of $\eta_1(a_1+b_1\omega)+\mu_1-\eta_2(a_2+b_2\omega)-\mu_2$ was equal to the identity polynomial, which takes forever. Sympy uses the first algorithm by default and alternatively something involving Gröbner bases which I did not understand.
Of course the object-oriented-programming approach to solving this problem is declaring a class which represents an algebraically closed structure that holds all the numbers you're gonna need (and you better foresee all of them!). For example in c++:
// a member of a real or imaginary quadratic field
// a+b*sqrt(n) where a,b,n are integers
class quadratic_int {
int a;
int b;
int n;
quadratic_int operator+(quadratic_int p) const {
// implement addition
}
// ...
};
This will be fast to perform operations on, but if you have a number which is not of the form $a+b\sqrt n$ you need to rewrite everything from scratch, plus if $n\equiv 1 \pmod 4$, $a$ and $b$ need to be half integers, and if you wanted to perform operations on the field $\mathbb Q[\sqrt 2]$ instead the class should probably use a template type instead of ints...all of this complexity means that by the end the class will be over a hundred lines long.
If we could compute minimal polynomials fast though, we could have a class called algebraic which internally represents a number as its minimal polynomial together with a way of distinguishing it from the other roots. Then its usage could be something like:
algebraic a {"x^2-2", 0}; // the 0th root of x^2-2 i.e. sqrt(2)
algebraic b {"x^3-3", 0}; // in 3^(1/3)
// x = 2^(1/2)+3^(1/3)
// so x is an algebraic number of degree 6
// the internal representation of x holds the polynomial
// t^6-6t^4-6t^3+12t^2-36t+1
// as a list of 7 integers (1,0,-6,-6,12,-36,1) plus the number of the root
// which is an integer from 0 to 5
auto x = a + b;
Where the 0th, 1st, etc roots are ordered say by absolute value and then by complex argument.
This way we still get fast arithmetic without having to rewrite our class every time we need to work on a different algebraic structure.
Statement of the question:
What is the fastest algorithm to calculate the minimal polynomial of a number expressed in radicals? It should have good asymptotic complexity in the number of terms/the degree of the input. Performance for small entries is not so important because a look-up table can be implemented.
