# Out of memory when computing a coefficient of a large symmetric polynomial

 f[b_, c_, d_] := (-b - c - d + 1)^1024*(b - c - d + 1)^1024*
(-b + c - d + 1)^1024*(b + c - d + 1)^1024*(-b - c + d +
1)^1024*
(b - c + d + 1)^1024*(-b + c + d + 1)^1024*(b + c + d +
1)^1024;
Coefficient[f[b, c, d], b^12*c^8*d^6]


The above code runs for a while, before the kernel crashes due to being out of memory. What can be done in order to complete the computation successfully?

The polynomial is indeed large, but it is also symmetric. Can the symmetry be exploited in order to complete the computation successfully?

• I get answer is -841854147465097693828672676093256048820224. It needed about 26 Gb of RAM to use and took about 3 minutes to run.. Feb 5, 2022 at 23:20
• I don't have enough RAM on my machine currently, but good to know that the difference can be overcome. Feb 5, 2022 at 23:25
• It crashed on my iMac (16Gb of RAM) but on my workstation (1Tb of RAM) obviously no problem. Feb 5, 2022 at 23:29
• Two ideas: Could you afford a small fast solid state hard drive, maybe 64 GB, and set your computer up to use that as 64 GB of swap space and try the calculation. That may be 100-1000 times slower than ram, but you can see if it works. OR could you rewrite the code to be an 8 level deep nested For loop where it would calculate each term of the polynomial, but not save all those, and only add to a running total what the final coefficient would be. That would again be far slower, but use far less memory if you coded it carefully. Mathematica sometimes works best with WAY WAY too much memory.
– Bill
Feb 5, 2022 at 23:39
• Reported (belatedly) as a bug. Jul 6, 2023 at 17:53

Split it into 3 steps to get result immediately.

Coefficient[Coefficient[Coefficient[f[b, c, d], b^12], c^8], d^6]

(*   -841854147465097693828672676093256048820224   *)

• (+1) Kinda makes Coefficient[f[b, c, d], b^12*c^8*d^6] seem like a bug. (Alternative coding: Fold[Coefficient, f[b, c, d], {b^12, c^8, d^6}]) Feb 6, 2022 at 16:22
• Was there some kind of general principle you followed to get this restatement of the problem or was it a brainstorm or something else? Feb 7, 2022 at 0:14
• @davidbak The consideration behind is simple: In the result of Coefficient[f[b, c, d], b^12]  all terms with c^8 are of course contained, the same with d^6. The next Coefficient comand selects it and so on. The only task was to test it then. Let me emphasize, this works so well, because ot the very effictive routines of Coefficient for one coefficient. Also see the comment of @MichaelE2, who regards failure of multiple selection as a bug. Feb 7, 2022 at 7:08

To get the coefficient of $$b^\beta c^\gamma d^\delta$$, for $$\beta=12$$, $$\gamma=8$$, $$\delta=6$$, we take advantage of the fact that the coefficients $$(poly)^n$$ are computed by the Multinomial[] theorem. To apply it, we need to Solve[] for the combinations of monomials in $$poly$$ whose product is the monomial $$b^\beta c^\gamma d^\delta$$.

base = Expand[poly[b, c, d] // #^(1/1024) & // PowerExpand];
monomialRules =   (* label the 35 monomials of $$poly$$ X[k] *)
MapIndexed[X@First@#2 &, List @@ base] -> List @@ base //
Thread; (* the 1/1024-th root yields $$poly$$ *)
Block[{β = 12, γ = 8, δ = 6}, (* pick powers *)
combos = Solve[ (* solve for the combinations of monomials *)
{Sum[q[k], {k, 35}] == 1024,
Sum[Exponent[X[k] /. monomialRules, b] q[k], {k, 35}] == β,
Sum[Exponent[X[k] /. monomialRules, c] q[k], {k, 35}] == γ,
Sum[Exponent[X[k] /. monomialRules, d] q[k], {k, 35}] == δ,
And @@ Thread[Array[q, 35] >= 0]},
Array[q, 35], Integers]
];

• When I evaluate your code I get a sum of the form Multinomial[q[1],...,q[35]]*q[1]*X[1]+...+Multinomial[q[1],...,q[35]]*q[35]*X[35], instead of the expected numerical output. I suppose that the q[i] and X[i] are not getting solved for some reason. Feb 6, 2022 at 8:58
• @PalmTopTigerMO I had a % from a missing line. Fixed now. Feb 6, 2022 at 15:45