11
$\begingroup$ 
 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?

Improve this question
$\endgroup$
5
  • 3
    $\begingroup$ I get answer is -841854147465097693828672676093256048820224. It needed about 26 Gb of RAM to use and took about 3 minutes to run.. $\endgroup$
    – 1729taxi
    Feb 5, 2022 at 23:20
  • $\begingroup$ I don't have enough RAM on my machine currently, but good to know that the difference can be overcome. $\endgroup$
    – PalmTopTigerMO
    Feb 5, 2022 at 23:25
  • 2
    $\begingroup$ It crashed on my iMac (16Gb of RAM) but on my workstation (1Tb of RAM) obviously no problem. $\endgroup$
    – 1729taxi
    Feb 5, 2022 at 23:29
  • 2
    $\begingroup$ 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. $\endgroup$
    – Bill
    Feb 5, 2022 at 23:39
  • 1
    $\begingroup$ Reported (belatedly) as a bug. $\endgroup$
    – Daniel Lichtblau
    Jul 6 at 17:53

2 Answers 2

Reset to default
19
$\begingroup$

Split it into 3 steps to get result immediately.

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

(*   -841854147465097693828672676093256048820224   *)
Improve this answer
$\endgroup$
3
  • 3
    $\begingroup$ (+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}]) $\endgroup$
    – Michael E2
    Feb 6, 2022 at 16:22
  • 1
    $\begingroup$ Was there some kind of general principle you followed to get this restatement of the problem or was it a brainstorm or something else? $\endgroup$
    – davidbak
    Feb 7, 2022 at 0:14
  • 1
    $\begingroup$ @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. $\endgroup$
    – Akku14
    Feb 7, 2022 at 7:08
10
$\begingroup$

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]
  ];
monomialCoeffs = Thread[
   Array[X, 35] ->
       (CoefficientList[Array[X, 35] /. monomialRules, {b, c, d}] // 
      Map@Flatten // Map@DeleteCases[0] // Flatten)
   ];
(Multinomial @@ Array[q, 35] /. combos) *
 (Array[X, 35]^Array[q, 35] /. monomialCoeffs /. combos // Map@Apply@Times) //
  Total

(*  -841854147465097693828672676093256048820224  *)
Improve this answer
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – PalmTopTigerMO
    Feb 6, 2022 at 8:58
  • 1
    $\begingroup$ @PalmTopTigerMO I had a % from a missing line. Fixed now. $\endgroup$
    – Michael E2
    Feb 6, 2022 at 15:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.