Finding coefficients in polynomial efficiently

Question

Given a transitive $G$-set $ M $. I'm interested in finding the number of fixed points of $ G $ acting on $\operatorname{Pot}_a( M ) := \left\{ N\subseteq M; |N| = a \right\} $ by using the table of marks of $G$. Let $U_1, \dots, U_s $ be representatives of the conjugacy classes of subgroups of $ G $. For each of them, a polynomial

$$p_{U_i} := \prod_{l=1}^s \left( \sum_{j=0}^{m_{i,l}} \binom{m_{i,l}}{j}\cdot x^{\lambda(i,l)\cdot j} \right)$$

is defined, where $\lambda(i,l), m_{i,l}\in \mathbb N_0$. Let $ p := \sum_{i=1}^s p_i $. The coefficient of $x^p$ is the number of fixed points $G$ has on $\operatorname{Pot}_p( M ) $. The problem is to calculate these coefficients efficiently. My approach was to calculate

$$p_{U_i} := \prod_{l=1}^s \left( \sum_{j=0}^{m_{i,l}} \binom{m_{i,l}}{j}\cdot x^{\lambda(i,l)\cdot j} \operatorname{mod}x^{p+1} \right) \operatorname{ mod } x^{p+1}.$$

Each transitive $G$-set corresponds to a row in the table of marks of $G$ since $M\cong G/U_k$ for some $U_k \leq G$. It turns out that calculating the fixedpoints of $S_8$ on $\operatorname{Pot}_2( G/U_p )$ for $p=2 \dots 148$ takes more than 10 hours ( still not finished ) which is not acceptable. To define the polynomial I used the internal GAP functions Product, Sum and mod. Do you have any suggestions how I can calculate the coefficients faster?

Answer 1

Thanks, that's really an interesting question!

$0$-th approach

When multiplying polynomials with integer coefficients, the absolute value of coefficients may increase quite fast. It's important to check that GAP is compiled with GMP library for fast integer arithmetic. You should see gmp in the 4-th line of the GAP startup screen:

 Libs used:  gmp, readline

GMP support was introduced in GAP 4.5 (2012) and is now a default option, so most likely this will be the case.

Step 1

The first step is to directly assemble the polynomial from the list of coefficients instead of calculating it as the sum of monomials. In the following example I create a random list l of a $1000$ integers between $-10$ and $10$, and then compare these two ways of creating a polynomial with coefficients given by this list:

gap> x:=Indeterminate(Rationals,"x");
gap> l:=List([1..10000],i->Random(Integers));;
gap> f:=Sum([1..10000], i-> l[i]*x^i);;time;
2268
gap> f1:=LaurentPolynomialByCoefficients(FamilyObj(1),l,1);;time;
1
gap> f=f1;
true

time returns the CPU time of the last command in milliseconds, so you see that LaurentPolynomialByCoefficients does this instantly, while the other approach takes more than two seconds. See more in "Creation of Rational Functions" from the chapter Polynomials and Rational Functions of the GAP reference manual - there are also more low-level functions like PolynomialByExtRep and PolynomialByExtRepNC, LaurentPolynomialByExtRep and LaurentPolynomialByExtRepNC, where NC stands for "no check" and only should be used if speed is required and the arguments are known to be in correct form. Unless this is guaranteed for the parameters, LaurentPolynomialByCoefficients should be used to be on the safe side.

Obviously, there is now no need in mod at this step any more, since one could simply truncate the list of coefficients (or not to compute the full list at all).

Remark: Step 1 has been verified (see question author's comments above) and reduced in one example the overall runtime from 22 hours to 59 minutes. Steps below are just suggestions, requiring more experimenting and fine-tuning the code.

Step 2(a)

Perhaps this could be improved further by using the functionality described in Arithmetic for External Representations of Polynomials or even Vectors as coefficients of polynomials. The main goal is to avoid back and forth conversions between polynomials as GAP objects and polynomials given by their external representations.

Added later: in particular, consider ProductCoeffs which allows you to multiply polynomials only partially: ProductCoeffs( list1[, len1], list2[, len2] ) assumes that $p1$ (and $p2$) are polynomials given by the first len1 (len2) entries of the coefficient list list2 (list2). If len1 and len2 are omitted, they default to the lengths of list1 and list2. Then it returns the coefficient list of the product of $p1$ and $p2$.

If lengths are not restricted, the performance is the same:

gap> x:=Indeterminate(Rationals,"x");;
gap> deg:=10000;;
gap> l1:=List([1 .. deg],i->Random(Integers));;
gap> f1:=LaurentPolynomialByCoefficients(FamilyObj(1),l1,0);;
gap> l2:=List([1 .. deg],i->Random(Integers));;
gap> f2:=LaurentPolynomialByCoefficients(FamilyObj(1),l2,0);;
gap> f:=f1*f2;;time;
3819
gap> t:=ProductCoeffs(CoefficientsOfUnivariatePolynomial(f1), CoefficientsOfUnivariatePolynomial(f2));;time;
3787
gap> CoefficientsOfUnivariatePolynomial(f)=t;
true

But now assume that we want to calculate coefficient for $x^{10}$ in the product of $p1$ and $p2$. Then we do not have to use monomials of degree 11 and higher. Thus we do:

gap> p:=10;
10
gap> tcut:=ProductCoeffs(CoefficientsOfUnivariatePolynomial(f1), p+1, CoefficientsOfUnivariatePolynomial(f2), p+1);;time;
0
gap> CoefficientsOfUnivariatePolynomial(f){[1..p+1]}=tcut{[1..p+1]};
true

so we are getting the result instantly, skipping dealing with all monomials of higher degrees. In the same setup, I've got $43$ms for $p=1000$ and $977$ms for $p=5000$, so your mileage may vary, and the smaller $p$ is, the better should be the speedup.

Step 2(b)

One could try to experiment with parallelising multiplication of polynomials (given by internal or external representations). The number of polynomials to multiply is the number of conjugacy classes of $G$, e.g. $22$ for $S_8$. The simplest approach would be to split the list of polynomials into chunks and then find a product of each chunk in parallel, and then multiply the resulting products (possibly with another parallel iteration, depending on the number of chunks and nodes). Anyhow, the degree will be rising, and the obvious bottleneck would be the last step with the sequential multiplication of two or three polynomials. If you wish to try to experiment with this, look at the SCSCP package. I could also provide a code for parallel Karatsuba multiplication of polynomials with SCSCP package - for large degree polynomials, it might also speed it up a bit.

Step 2(c)

Another speculative suggestion comes from the fact that you may be interested only in a single coefficient. That would be much more helpful if polynomials would be sparse, but still it may be worth to try in the dense case. The idea would be to iterate over the cartesian product of degrees of monomials (using IteratorOfCartesianProduct), check for each tuple whether the sum of degrees of monomials is equal to the degree in which you're interested in, and if so then add the product of the respective coefficients to some accumulating variable. The pseudo-code to calculate p-th coefficient in GAP could look like

degs:= < list of degrees of polynomials >
degranges := List( degs, i -> [1..i] ); # list of ranges 
# for sparse polynomial, take the list of degrees of monomials instead of a range

coeff := 0;

iter:=IteratorOfCartesianProduct(degranges);

for t in iter do
  # t is a tuple of degrees of monomials 
  if Sum(t) = p then
    coeff := coeff + Product( List( [1..Length(t)], i -> 
                              < coefficient for x^t[i] in the i-th polynomial > 
  fi;
od;