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][1]

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

![enter image description here][2]

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},
      
      (* fugly patch for now *)
      If[(Tr[ipow*epow] < Tr@target) || (Min@target < Min@ipow), Return[0]];

      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]}, 
            (* Another fugly patch for non-conformance *)
             If[Length@tmp < left, Return[{{}}]];
             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](https://mathematica.stackexchange.com/questions/46261/undocumented-form-for-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).

Update 25/02/2015 per comments:

Looking at your test harness, there's an obvious issue - you can't just *Quiet* something and expect that errors sans messages means the same performance. Even though using *Quiet* does usually reduce timings when errors cause messages, there's still overhead in *Mathematica* handling the issue, and of course the likelihood that there's a cascade of problems presented to the code... also on large result lists, using rules can be quite inefficient for replacements. Look into *Dispatch* and using *Replace* directly with a limited level specification if you have tasks that need to do such things in the future.

Since the original question, comments, *Math.SE* question, etc. all implied the feeding of well-conditioned exponent lists, this was built with that in mind. I've added another ugly patch for now, but be advised, it just short-ciruits the problem, but it also means there's work, sometimes a lot, getting done where it need not be.

That said, I took your code, and just changed one line to be:

    getCoeff[Sequence @@ Transpose@Tally@list1, #] & /@ list2

in your *partialONk* function. No other changes made, so there's extra work converting the argument format, and the extra work being done in my function because bad sets are passed to it. Nonetheless, the following is the result (netbook timings, I ran out of patience at 12):

![enter image description here][3]

Both grow quickly (expected from the problem), but the original *coeff* *explodes*.

There appears to be *much* overhead from the rest of your code, based on just timing the coefficient functions. I'd venture re-thinking it, and in particular generating only valid exponent sets, will go a long way in improving timing. As it stands, much of the boost in finding the coefficients is getting swamped in other overhead. By example, the two functions, over a range of 1 to 5 components, each with random powers between 2 and 5 for internal/external power, and a random exponent set of appropriate size. I gave up at six components - *getCoeff* finished in under a second, I aborted *coeff* after twenty minutes. At five components, *getCoeff* was ~2800 times faster... N.B. This had to be done with a *Log scale* so *getCoeff* would even show up on the plot:

![enter image description here][4]


  [1]: https://i.sstatic.net/7HHaq.png
  [2]: https://i.sstatic.net/FLzex.png
  [3]: https://i.sstatic.net/XT4eO.png
  [4]: https://i.sstatic.net/I8PUF.png