#include <utility>
#include <vector>
std::vector< std::pair<int, int> > factor_table;
void fill_sieve( int n )
{
factor_table.resize(n+1);
for( int i = 1; i <= n; ++i )
factor_table[i] = std::pair<int, int>(i, 1);
for( int j = 2, j2 = 4; j2 <= n; (j2 += j), (j2 += ++j) ) {
if (factor_table[j].second == 1) {
int i = j;
int ij = j2;
while (ij <= n) {
factor_table[ij] = std::pair<int, int>(j, i);
++i;
ij += j;
}
}
}
}
std::vector<unsigned> powers;
template<int dir>
void factor( int num )
{
while (num != 1) {
powers[factor_table[num].first] += dir;
num = factor_table[num].second;
}
}
template<unsigned N>
void calc_combinations(unsigned (&bin_sizes)[N])
{
using std::swap;
powers.resize(0);
if (N < 2) return;
unsigned& largest = bin_sizes[0];
size_t sum = largest;
for( int bin = 1; bin < N; ++bin ) {
unsigned& this_bin = bin_sizes[bin];
sum += this_bin;
if (this_bin > largest) swap(this_bin, largest);
}
fill_sieve(sum);
powers.resize(sum+1);
for( unsigned i = largest + 1; i <= sum; ++i ) factor<+1>(i);
for( unsigned bin = 1; bin < N; ++bin )
for( unsigned j = 2; j <= bin_sizes[bin]; ++j ) factor<-1>(j);
}
#include <iostream>
#include <cmath>
int main(void)
{
unsigned bin_sizes[] = { 5, 15 };
calc_combinations(bin_sizes);
char* sep = "";
for( unsigned i = 0; i < powers.size(); ++i ) {
if (powers[i]) {
std::cout << sep << i;
sep = " * ";
if (powers[i] > 1)
std::cout << "**" << powers[i];
}
}
std::cout << "\n\n";
}
I2luY2x1ZGUgPHV0aWxpdHk+CiNpbmNsdWRlIDx2ZWN0b3I+CgpzdGQ6OnZlY3Rvcjwgc3RkOjpwYWlyPGludCwgaW50PiA+IGZhY3Rvcl90YWJsZTsKdm9pZCBmaWxsX3NpZXZlKCBpbnQgbiApCnsKCWZhY3Rvcl90YWJsZS5yZXNpemUobisxKTsKCWZvciggaW50IGkgPSAxOyBpIDw9IG47ICsraSApCgkJZmFjdG9yX3RhYmxlW2ldID0gc3RkOjpwYWlyPGludCwgaW50PihpLCAxKTsKCWZvciggaW50IGogPSAyLCBqMiA9IDQ7IGoyIDw9IG47IChqMiArPSBqKSwgKGoyICs9ICsraikgKSB7CgkJaWYgKGZhY3Rvcl90YWJsZVtqXS5zZWNvbmQgPT0gMSkgewoJCQlpbnQgaSA9IGo7CgkJCWludCBpaiA9IGoyOwoJCQl3aGlsZSAoaWogPD0gbikgewoJCQkJZmFjdG9yX3RhYmxlW2lqXSA9IHN0ZDo6cGFpcjxpbnQsIGludD4oaiwgaSk7CgkJCQkrK2k7CgkJCQlpaiArPSBqOwoJCQl9CgkJfQoJfQp9CgpzdGQ6OnZlY3Rvcjx1bnNpZ25lZD4gcG93ZXJzOwoKdGVtcGxhdGU8aW50IGRpcj4Kdm9pZCBmYWN0b3IoIGludCBudW0gKQp7Cgl3aGlsZSAobnVtICE9IDEpIHsKCQlwb3dlcnNbZmFjdG9yX3RhYmxlW251bV0uZmlyc3RdICs9IGRpcjsKCQludW0gPSBmYWN0b3JfdGFibGVbbnVtXS5zZWNvbmQ7Cgl9Cn0KCnRlbXBsYXRlPHVuc2lnbmVkIE4+CnZvaWQgY2FsY19jb21iaW5hdGlvbnModW5zaWduZWQgKCZiaW5fc2l6ZXMpW05dKQp7Cgl1c2luZyBzdGQ6OnN3YXA7CgoJcG93ZXJzLnJlc2l6ZSgwKTsKCWlmIChOIDwgMikgcmV0dXJuOwoKCXVuc2lnbmVkJiBsYXJnZXN0ID0gYmluX3NpemVzWzBdOwoJc2l6ZV90IHN1bSA9IGxhcmdlc3Q7Cglmb3IoIGludCBiaW4gPSAxOyBiaW4gPCBOOyArK2JpbiApIHsKCQl1bnNpZ25lZCYgdGhpc19iaW4gPSBiaW5fc2l6ZXNbYmluXTsKCQlzdW0gKz0gdGhpc19iaW47CgkJaWYgKHRoaXNfYmluID4gbGFyZ2VzdCkgc3dhcCh0aGlzX2JpbiwgbGFyZ2VzdCk7Cgl9CglmaWxsX3NpZXZlKHN1bSk7CgoJcG93ZXJzLnJlc2l6ZShzdW0rMSk7Cglmb3IoIHVuc2lnbmVkIGkgPSBsYXJnZXN0ICsgMTsgaSA8PSBzdW07ICsraSApIGZhY3RvcjwrMT4oaSk7Cglmb3IoIHVuc2lnbmVkIGJpbiA9IDE7IGJpbiA8IE47ICsrYmluICkKCQlmb3IoIHVuc2lnbmVkIGogPSAyOyBqIDw9IGJpbl9zaXplc1tiaW5dOyArK2ogKSBmYWN0b3I8LTE+KGopOwp9CgojaW5jbHVkZSA8aW9zdHJlYW0+CiNpbmNsdWRlIDxjbWF0aD4KaW50IG1haW4odm9pZCkKewoJdW5zaWduZWQgYmluX3NpemVzW10gPSB7IDUsIDE1IH07CgljYWxjX2NvbWJpbmF0aW9ucyhiaW5fc2l6ZXMpOwoJY2hhciogc2VwID0gIiI7Cglmb3IoIHVuc2lnbmVkIGkgPSAwOyBpIDwgcG93ZXJzLnNpemUoKTsgKytpICkgewoJCWlmIChwb3dlcnNbaV0pIHsKCQkJc3RkOjpjb3V0IDw8IHNlcCA8PCBpOwoJCQlzZXAgPSAiICogIjsKCQkJaWYgKHBvd2Vyc1tpXSA+IDEpCgkJCQlzdGQ6OmNvdXQgPDwgIioqIiA8PCBwb3dlcnNbaV07CgkJfQoJfQoJc3RkOjpjb3V0IDw8ICJcblxuIjsKfQo=