Skip to main content
2 of 5
Added finder code
ciao
  • 26k
  • 2
  • 61
  • 142

This is not (at this time) a direct answer to the core question, perhaps better thought of as a placeholder. Nonetheless, it will serve as an example of how to speed up this problem, and like things in general.

You need to start thinking in "Mathematica" terms. That is, whenever possible, think of how you might manipulate things en masse, as operations on vectors, matrices, lists, etc. If you can vectorize things in Mathematica, it;s nearly always a win.

By way of example, I started to decode your code in the OP to get a better understanding of what you're after. Here's the first function, Mathematica-ized:

newip[n_] := Times @@@ (Prime@Range[Length /@ #]^#) &[IntegerPartitions[n]];

How much of a difference does that make?

enter image description here

We have to view it in Log scale to even see the faster way:

enter image description here

By the time N is 30, this is over 250X faster.

I'll update this as I go through your code, and if the current ideas for the core problem flesh out, I'll add them.

Update 24/02/2015:

OK, here's the coefficient finder I've been tooling with.

Some caveats: It's obviously a work-in-progress, parts within that could be combined are separated to allow profiling of them individually, some parts certainly could be further optimized (e.g. using outer vs mapping over sets for low cardinality partition sets), etc., but perhaps there's some useful stuff here for you.

I've not tested it on huge arguments compared to native Mathematica methods - I do my play puzzling on leisure netbook & laptop, so no testing on real workstations was done...

This takes arguments similar to the earlier answer in the linked question, e.g.

(x[1] + x[2] + x[3] + x[4] + x[5] + x[6])^10 (x[1]^2 + x[2]^2 + 
   x[3]^2 + x[4]^2 + x[5]^2 + x[6]^2)^10 (x[1]^3 + x[2]^3 + x[3]^3 + 
   x[4]^3 + x[5]^3 + x[6]^3)^2

is represented with the first two arguments being the "internal" and "external" powers:

{1, 2, 3, 2, 2} , {10, 1, 2, 4, 5}

the third argument is the exponent set you're after, e.g.:

{17, 9, 7, 1, 1, 1}

representing an exponent set in six variables. You can use things like

{17, 9, 10} or  {17, 9, 10, 0, 0, 0} 

for terms with zero exponents, they are treated equivalently.

A quick test on the above case gives this an ~2500X speed advantage over your coeff, and about a 10X advantage over Coefficient in Mathematica. I'd expect the gap to grow as problem size grows. It is also much gentler on memory resources.

As stated in my comment, getting all is still faster using something like CoefficientArray (I've yet to find the Kryptonite for this problem, any luck for you over at Math.StackExchange?), but I'm guessing that's a no-go for you due to memory or other resource pressure?

getCoeff[ipow_, epow_, expset_] :=
 Module[{f, fastRF, myIP, px = {{}}, target = DeleteCases[expset, 0], 
   set = Join @@ MapThread[ConstantArray, {ipow, epow}], nres, x, bin,
    tset, breaks, accs, splits, seqs, res},
  
  fastRF[a_List, b_List] := Module[{c, o, x}, c = Join[b, a];
    o = Ordering[c];
    x = 1 - 2 UnitStep[-1 - Length[b] + o];
    x = FoldList[Max[#, 0] + #2 &, x];
    x[[o]] = x;
    Pick[c, x, -1]];
  
  myIP[set_, used_, left_, n_] := 
   With[{tmp = fastRF[set, used]}, 
         IntegerPartitions[n, Length@tmp - left, DeleteDuplicates@tmp]];
  
  tset = Tally[set];
  
  Module[{cnt = #[[2]], tmp, reset = #[[2]]}, 
     bin[#[[1]], x_] := (tmp = Binomial[cnt, x]; cnt -= x; 
       If[cnt == 0, cnt = reset]; tmp)] & /@ tset;
  
  (f[#[[1]], z_] = z <= #[[2]]) & /@ Tally[set];
  
  nres = Module[{z = 0}, 
    Nest[With[{zz = #}, z++; 
       Join @@ (Module[{mip = myIP[set, #, Length@target - z, target[[z]]], tmp = #},
                       mip = Join[tmp, #] & /@ mip;
            
            Pick[mip, Apply[And, Apply[f, Tally /@ mip, {2}], {1}]]] & /@ 
                             zz)] &, px, Length@target]];
  
  breaks = Accumulate@target;
  accs = Accumulate /@ nres;
  
  Reap[Do[
     x = Join @@ Position[accs[[z]], Alternatives @@ breaks];
     x = Differences[Prepend[x, 0]];
     splits = 
      Transpose@{Span @@@ Transpose[{Accumulate[Prepend[Most@x, 0]] + 1, Accumulate[x]}]};
     segs = Tally /@ Rest@Extract[nres[[z]], Prepend[splits, {{}}]];
     res = Times @@ Join @@ Apply[bin, segs, {2}];
     Sow[res];, {z, Length@nres}]] // 
   If[#[[2]] == {}, 0, Tr@#[[2, 1]]] &
  
  ]

A quick run-through:

Sub-function myIP is a restricted partition generator. It takes size and member restrictions and is used to build the candidates for each of the target exponent set components. It uses Mr. Wizard's fastRF to cull the restrictions. I may have a faster method, but it's an inconsequential part of the process, and his code is so pretty, so I kept it.

I the build a set of custom binomial coefficient functions bin[...]. These return the coefficient given the target value but keep track internally of how many of that value are available, resetting to the initial tally when exhausted for each candidate partition set.

Following, I build a set of f[...] functions that allow quick validation of the tallies of candidate partition sets (i.e., does the set have a valid amount of each possible piece).

The next step fills nres with valid partition sets using the above.

Finally, I parse the candidate partition sets, splitting them at the points where the running "stuttering" totals equal the target coefficients, tally these, extract the results using an undocumented form of extract I discovered for speed, and applying the adaptive bin binomial coefficient extractors. Lastly, I check if the result is empty (coefficient set was invalid) and return 0 (invalid) or the total of the coefficients (valid).

I hope at least some piece of this finds use, and will continue pondering this interesting puzzle (I'm convinced there's a better way, mathematically justified, that perhaps someone with deep multinomial/poly-product knowledge can answer with over at your Math SE query).

ciao
  • 26k
  • 2
  • 61
  • 142