polynomial multiplication is unexpectedly slow
Hello!
I'm writing some code to solve by radicals a polynomial with solvable galois group, and an important part of doing this is the special case of a cyclic field extension.
I've implemented the algorithm (probably naively) and while it _works_, it's extremely slow to do a polynomial multiplication in the middle part of the algorithm. Granted, we end up with a polynomial of degree $n!$, but for $n=5$ I feel like it shouldn't be quite as slow as it is (though maybe I'm wrong). As an aside, if anybody knows of a more efficient algorithm, I would love to hear about it. I'm thinking of looking through the references of this paper, but I haven't had time yet.
Here's a minimal excerpt from the code which shows the part that runs slowly:
f = x^3 + x^2 - 2*x - 1
n = f.degree(x)
# make a number field where w is a nth root of 1
K.<w> = CyclotomicField(n)
# make vars for the roots.
# a[i] = alpha_i
R = PolynomialRing(K, n, 'a')
a = R.gens()
v = sum([w^k * a[k] for k in range(n)])^3
# build the conjugates of v.
# it turns out we only need one from each conjugacy class,
# which keeps the degree down.
Sn = SymmetricGroup(n)
conjugates = []
for tau in Sn:
if v * tau not in conjugates:
conjugates += [v * tau]
S.<Y> = PolynomialRing(R)
psi = 1
for c in conjugates:
psi *= (Y - c)
It took my desktop (Ryzen 5 3600) running overnight to fully compute this polynomial. I've tried using prod(Y-c for c in conjugates) as well as product, rather than the explicit loop I have here. I've also tried making Y a symbolic variable rather than an element of a polynomial ring, but nothing has appreciably decreased the runtime.
The first 10 multiplications or so happen without much hassle, but very quickly things start going very slowly. If anyone has any ideas, I would love to hear them!
Thanks in advance ^_^
Probably you should use somewhere that the
a[i]are roots; currently they don't satisfy any non-trivial polynomial equations. Also you can work in a polynomial ring over Q if you add an extra variable for the root of unity. In any case, I would expect to seereducesomewhere.@rburing -- that's exactly the next step. We know that the coefficients of
fare symmetric polynomials in the rootsa[i]. The problem is thatv_isn't_ symmetric. However, once we multiply by all the conjugates of the symmetric group, then the coefficients ofpsi_are_ symmetric, and we can write them in terms of the elementary symmetric polynomials (that is, the coefficients off). The references I've seen have all needed this step before we can actually evaluate the coefficients ofpsidown to constants, but there might be a way to avoid it? As for working in a polynomial ring over Q, is that more efficient? I'm thinking now about evaluating each coefficient by hand and reducing them one by one, rather than doing a full polynomial multiply and reducing after.It's probably a good exercise to write down the expected number of terms in the coefficients of this degree 120 polynomial. Sure, a degree 120 polynomial over QQ usually quite workable (if the height of the coefficients remains moderate), but in your case the coefficients are polynomials themselves. If you multiply such polynomials together, you can get very many terms very quickly, and with a lot of cross terms. A few things you could try: - use a balanced multiplication technique (try to multiply low-degree polynomials with low-degree ones, so that there's only one very expensive multiplication) - flatten your polynomial ring to hopefully optimize the polynomial multiplications a bit. I suspect you're just using a very computationally intensive method.
Sis a nested polynomial ring, for which multiplication is very slow in Sage. It uses the Karatsuba algorithm, which I think is a particularly bad choice if the coefficients are polynomials themselves.Instead, you can compute the product in
S.flattening_morphism().codomain(), the flattened polynomial ring @nbruin suggested, and then convert the result back toSusingS.flattening_morphism().section(). This might be considerably faster, since it uses libsingular.Thank you all for your ideas! I've modified my code to use the
S.flattening_morphism().codomain(), and this has helped significantly (though not enough to get the degree $5$ case to terminate). As for using a balanced multiplication technique, I thinkprodalready does that. According toprod??:prod = balanced_list_prod(x, 0, n, recursion_cutoff). I suppose I can try to implement my own balanced product, but I feel like anything I do is unlikely to be better than the official version. Do you think it would be helpful if I edited my question to show the entire script I'm trying to run, rather than just this minimal part of the example? Computingpsiis still the slowest part of the process, but perhaps I can compute it in stages, etc. based on what I'm doing next? Thanks again!@mwageringel, @nbruin, @rburing
As an aside, I wasn't able to get this particular implementation working efficiently, but I found an algorithm to more efficiently compute the object I was trying to construct in this question. In case it helps anyone else, you can find the details on my blog.