Problem 475: Music Festival
View on Project EulerProject Euler Problem 475 Solution
EulerSolve provides an optimized solution for Project Euler Problem 475, Music Festival, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Let \(12n\) musicians be split on the first day into \(m=3n\) quartets. On the second day they must be split into \(t=4n\) trios, and no pair of musicians is allowed to play together twice. Because every pair inside a first-day quartet has already met, a valid second-day trio can contain at most one musician from each quartet. The task is to count all valid second-day schedules modulo $$10^9+7.$$ Mathematical Approach The fixed first-day partition gives \(m\) distinguishable quartets. The whole problem becomes: how many ways can \(t\) trios be formed so that every trio draws one musician from three different quartets, while each quartet contributes all four of its musicians exactly once? Step 1: Encode the schedule by an incidence matrix Introduce an \(m\times t\) matrix \(A\) with entries in \(\{0,1\}\): $$A_{q,r}=1 \iff \text{trio }r\text{ contains one musician from quartet }q.$$ The no-repeat condition forces \(A_{q,r}\) to be binary, because a trio cannot take two musicians from the same quartet. Each trio has exactly three musicians, so every column sum is $$\sum_{q=1}^{m} A_{q,r}=3.$$ Each quartet contains four musicians and each of them must appear on day 2 exactly once, so every row sum is $$\sum_{r=1}^{t} A_{q,r}=4.$$ Therefore the first combinatorial subproblem is to count \(0\)-\(1\) matrices with row sum \(4\) and column sum \(3\)....
Detailed mathematical approach
Problem Summary
Let \(12n\) musicians be split on the first day into \(m=3n\) quartets. On the second day they must be split into \(t=4n\) trios, and no pair of musicians is allowed to play together twice. Because every pair inside a first-day quartet has already met, a valid second-day trio can contain at most one musician from each quartet. The task is to count all valid second-day schedules modulo
$$10^9+7.$$
Mathematical Approach
The fixed first-day partition gives \(m\) distinguishable quartets. The whole problem becomes: how many ways can \(t\) trios be formed so that every trio draws one musician from three different quartets, while each quartet contributes all four of its musicians exactly once?
Step 1: Encode the schedule by an incidence matrix
Introduce an \(m\times t\) matrix \(A\) with entries in \(\{0,1\}\):
$$A_{q,r}=1 \iff \text{trio }r\text{ contains one musician from quartet }q.$$
The no-repeat condition forces \(A_{q,r}\) to be binary, because a trio cannot take two musicians from the same quartet.
Each trio has exactly three musicians, so every column sum is
$$\sum_{q=1}^{m} A_{q,r}=3.$$
Each quartet contains four musicians and each of them must appear on day 2 exactly once, so every row sum is
$$\sum_{r=1}^{t} A_{q,r}=4.$$
Therefore the first combinatorial subproblem is to count \(0\)-\(1\) matrices with row sum \(4\) and column sum \(3\).
Step 2: Track only how many times each quartet has already been used
Process the trio columns from left to right. After some number of processed trios, every quartet has been used \(0\), \(1\), \(2\), \(3\), or \(4\) times. Quartets already used \(4\) times are finished and never matter again, so the state only needs
$$\mathbf{c}=(c_0,c_1,c_2,c_3),$$
where \(c_j\) is the number of quartets used exactly \(j\) times so far. The number used \(4\) times is implicit:
$$c_4=m-c_0-c_1-c_2-c_3.$$
Initially no day-2 trio has been built, hence
$$\mathbf{c}_{\text{start}}=(m,0,0,0).$$
Step 3: Count one transition
To build the next trio, choose
$$x_j \text{ quartets from class }j \qquad (j=0,1,2,3),$$
subject to
$$x_0+x_1+x_2+x_3=3,\qquad 0\le x_j\le c_j.$$
This means the trio takes one musician from \(x_0\) previously unused quartets, one musician from \(x_1\) quartets already used once, and so on. The number of ways to make that choice is
$$\binom{c_0}{x_0}\binom{c_1}{x_1}\binom{c_2}{x_2}\binom{c_3}{x_3}.$$
There are only
$$\binom{3+4-1}{4-1}=\binom{6}{3}=20$$
possible quadruples \((x_0,x_1,x_2,x_3)\), so each state has a very small transition menu.
Step 4: Update the state and form the recurrence
Every chosen quartet moves from usage class \(j\) to usage class \(j+1\). Therefore
$$c_0'=c_0-x_0,$$
$$c_1'=c_1-x_1+x_0,$$
$$c_2'=c_2-x_2+x_1,$$
$$c_3'=c_3-x_3+x_2.$$
If \(D_s(c_0,c_1,c_2,c_3)\) denotes the number of ordered partial incidence matrices after \(s\) processed trios, then the recurrence is
$$D_{s+1}(c_0',c_1',c_2',c_3') \mathrel{+}= D_s(c_0,c_1,c_2,c_3)\prod_{j=0}^{3}\binom{c_j}{x_j}.$$
The base case is
$$D_0(m,0,0,0)=1.$$
After all \(t\) trio columns have been processed, every quartet must have been used four times, so the terminal state is
$$D_t(0,0,0,0).$$
This quantity is the number of ordered incidence matrices. Call it \(N_{\text{mat}}\).
Step 5: Convert matrix patterns into actual musician schedules
Once an incidence matrix is fixed, each quartet is attached to exactly four trio columns. Its four labeled musicians can be assigned to those four appearances in
$$4!=24$$
ways. Different quartets act independently, so the total refinement factor is
$$24^m.$$
However, the dynamic program processed the trio columns in an artificial order. A finished day-2 schedule is an unordered collection of \(t\) distinct trios, and every such schedule is counted once for each permutation of those \(t\) trio columns. Therefore we divide by \(t!\):
$$f(12n)=N_{\text{mat}}\cdot 24^m\cdot \frac{1}{t!}\pmod{10^9+7}.$$
Because \(10^9+7\) is prime, \((t!)^{-1}\) is computed with modular exponentiation using Fermat's little theorem.
Worked Example: \(12\) musicians
For \(12\) musicians we have \(n=1\), so
$$m=3,\qquad t=4.$$
There are exactly three first-day quartets. Every second-day trio must use three different quartets, so each trio must take one musician from each of the three quartets. Hence the incidence matrix is forced: every one of its \(3\times 4\) entries is \(1\). Thus
$$N_{\text{mat}}=1.$$
Now each quartet distributes its four musicians across the four trio positions in \(4!\) ways, giving
$$24^3$$
ordered assignments. Since the four trios themselves are unordered, divide by
$$4!=24.$$
So the final count is
$$f(12)=\frac{24^3}{24}=24^2=576,$$
which matches the checkpoint used by the implementation.
How the Code Works
The C++, Python, and Java implementations all follow the same plan. They compute \(m=3n\) and \(t=4n\), precompute the binomial values \(\binom{i}{0},\binom{i}{1},\binom{i}{2},\binom{i}{3}\) for \(0\le i\le m\), and enumerate the \(20\) admissible pick patterns \((x_0,x_1,x_2,x_3)\).
The dynamic program stores only the currently reachable states \((c_0,c_1,c_2,c_3)\) in a hash map, which keeps the state space compressed. For each processed trio, the implementation loops over the current states, applies every legal pick pattern, multiplies by the corresponding product of binomial coefficients, and accumulates the result for the next layer modulo \(10^9+7\).
After \(t\) layers, the implementation reads the value at the final state \((0,0,0,0)\), multiplies by \(24^m\), computes \(t!\) modulo \(10^9+7\), inverts that factorial by fast modular exponentiation, and multiplies once more. The small checkpoints \(f(12)=576\) and \(f(24)=509089824\) are used to verify that the recurrence has been implemented correctly.
Complexity Analysis
Let \(S\) be the maximum number of reachable DP states in any layer. Each state tries only \(20\) transition patterns, so the main loop runs in
$$O(20tS)=O(tS)$$
time. The binomial precomputation costs \(O(m)\), and the modular exponentiations for \(24^m\) and \((t!)^{-1}\) are negligible compared with the state transitions. The memory usage is
$$O(S),$$
because only the current and next hash-map layers are stored.
Footnotes and References
- Problem page: https://projecteuler.net/problem=475
- Incidence matrix: Wikipedia — Incidence matrix
- Dynamic programming: Wikipedia — Dynamic programming
- Binomial coefficient: Wikipedia — Binomial coefficient
- Fermat's little theorem: Wikipedia — Fermat's little theorem
Problem 475 source code
C++
#include <array>
#include <cstdint>
#include <iostream>
#include <string>
#include <unordered_map>
#include <vector>
namespace {
using u64 = std::uint64_t;
using u32 = std::uint32_t;
using i64 = std::int64_t;
constexpr u64 kMod = 1'000'000'007ULL;
struct Options {
int musicians = 600; // 12n
bool run_checkpoints = true;
};
bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& out) {
if (arg.rfind(prefix, 0U) != 0U) {
return false;
}
const std::string tail = arg.substr(prefix.size());
if (tail.empty()) {
return false;
}
try {
out = std::stoi(tail);
} catch (...) {
return false;
}
return true;
}
bool parse_arguments(int argc, char** argv, Options& options) {
for (int i = 1; i < argc; ++i) {
const std::string arg(argv[i]);
if (arg == "--skip-checkpoints") {
options.run_checkpoints = false;
continue;
}
if (parse_int_after_prefix(arg, "--musicians=", options.musicians)) {
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
if (options.musicians <= 0 || options.musicians % 12 != 0) {
std::cerr << "--musicians must be a positive multiple of 12.\n";
return false;
}
return true;
}
u64 mod_pow(u64 base, i64 exp) {
u64 out = 1ULL;
base %= kMod;
while (exp > 0) {
if (exp & 1LL) {
out = static_cast<u64>((__uint128_t)out * base % kMod);
}
base = static_cast<u64>((__uint128_t)base * base % kMod);
exp >>= 1LL;
}
return out;
}
u32 pack_state(const int c0, const int c1, const int c2, const int c3) {
return static_cast<u32>((c0 << 24) | (c1 << 16) | (c2 << 8) | c3);
}
void unpack_state(const u32 key, int& c0, int& c1, int& c2, int& c3) {
c0 = static_cast<int>((key >> 24) & 0xFFU);
c1 = static_cast<int>((key >> 16) & 0xFFU);
c2 = static_cast<int>((key >> 8) & 0xFFU);
c3 = static_cast<int>(key & 0xFFU);
}
u64 solve(const int musicians) {
const int n = musicians / 12;
const int m = 3 * n; // quartets on day 1
const int t = 4 * n; // trios on day 2
std::vector<std::array<u64, 4>> comb(static_cast<std::size_t>(m + 1));
for (int i = 0; i <= m; ++i) {
comb[static_cast<std::size_t>(i)][0] = 1ULL;
if (i >= 1) {
comb[static_cast<std::size_t>(i)][1] = static_cast<u64>(i);
}
if (i >= 2) {
comb[static_cast<std::size_t>(i)][2] = static_cast<u64>(i) * (i - 1) / 2ULL;
}
if (i >= 3) {
comb[static_cast<std::size_t>(i)][3] =
static_cast<u64>(i) * (i - 1) * (i - 2) / 6ULL;
}
}
std::vector<std::array<int, 4>> picks;
picks.reserve(20);
for (int x0 = 0; x0 <= 3; ++x0) {
for (int x1 = 0; x1 + x0 <= 3; ++x1) {
for (int x2 = 0; x2 + x1 + x0 <= 3; ++x2) {
const int x3 = 3 - x0 - x1 - x2;
picks.push_back({x0, x1, x2, x3});
}
}
}
std::unordered_map<u32, u64> cur;
std::unordered_map<u32, u64> nxt;
cur.reserve(200'000U);
nxt.reserve(200'000U);
cur[pack_state(m, 0, 0, 0)] = 1ULL;
for (int step = 0; step < t; ++step) {
nxt.clear();
for (const auto& kv : cur) {
int c0 = 0;
int c1 = 0;
int c2 = 0;
int c3 = 0;
unpack_state(kv.first, c0, c1, c2, c3);
const u64 ways_so_far = kv.second;
for (const auto& pick : picks) {
const int x0 = pick[0];
const int x1 = pick[1];
const int x2 = pick[2];
const int x3 = pick[3];
if (x0 > c0 || x1 > c1 || x2 > c2 || x3 > c3) {
continue;
}
const int nc0 = c0 - x0;
const int nc1 = c1 - x1 + x0;
const int nc2 = c2 - x2 + x1;
const int nc3 = c3 - x3 + x2;
u64 ways_pick = comb[static_cast<std::size_t>(c0)][static_cast<std::size_t>(x0)];
ways_pick = static_cast<u64>((__uint128_t)ways_pick *
comb[static_cast<std::size_t>(c1)]
[static_cast<std::size_t>(x1)] %
kMod);
ways_pick = static_cast<u64>((__uint128_t)ways_pick *
comb[static_cast<std::size_t>(c2)]
[static_cast<std::size_t>(x2)] %
kMod);
ways_pick = static_cast<u64>((__uint128_t)ways_pick *
comb[static_cast<std::size_t>(c3)]
[static_cast<std::size_t>(x3)] %
kMod);
const u64 add =
static_cast<u64>((__uint128_t)ways_so_far * ways_pick % kMod);
const u32 nkey = pack_state(nc0, nc1, nc2, nc3);
auto it = nxt.find(nkey);
if (it == nxt.end()) {
nxt.emplace(nkey, add);
} else {
u64 v = it->second + add;
if (v >= kMod) {
v -= kMod;
}
it->second = v;
}
}
}
cur.swap(nxt);
}
const u32 final_key = pack_state(0, 0, 0, 0);
const u64 matrix_count = cur[final_key];
const u64 assign_quartet_members = mod_pow(24ULL, m);
u64 fact_t = 1ULL;
for (int i = 2; i <= t; ++i) {
fact_t = static_cast<u64>((__uint128_t)fact_t * static_cast<u64>(i) % kMod);
}
const u64 inv_fact_t = mod_pow(fact_t, static_cast<i64>(kMod - 2ULL));
u64 out = static_cast<u64>((__uint128_t)matrix_count * assign_quartet_members % kMod);
out = static_cast<u64>((__uint128_t)out * inv_fact_t % kMod);
return out;
}
bool run_checkpoints() {
if (solve(12) != 576ULL) {
std::cerr << "Checkpoint failed: f(12)\n";
return false;
}
if (solve(24) != 509'089'824ULL) {
std::cerr << "Checkpoint failed: f(24)\n";
return false;
}
return true;
}
} // namespace
int main(int argc, char** argv) {
Options options;
if (!parse_arguments(argc, argv, options)) {
return 1;
}
if (options.run_checkpoints && !run_checkpoints()) {
return 1;
}
std::cout << solve(options.musicians) << '\n';
return 0;
}
Python
def solve():
MOD = 1000000007
musicians = 600
def mod_pow(base, exp):
r = 1; base %= MOD
while exp > 0:
if exp & 1: r = r * base % MOD
base = base * base % MOD; exp >>= 1
return r
n = musicians // 12
m = 3 * n # quartets
t = 4 * n # trios
comb = [[0]*4 for _ in range(m+1)]
for i in range(m+1):
comb[i][0] = 1
if i >= 1: comb[i][1] = i
if i >= 2: comb[i][2] = i*(i-1)//2
if i >= 3: comb[i][3] = i*(i-1)*(i-2)//6
picks = []
for x0 in range(4):
for x1 in range(4-x0):
for x2 in range(4-x0-x1):
x3 = 3-x0-x1-x2
picks.append((x0, x1, x2, x3))
def pack(c0,c1,c2,c3): return (c0<<24)|(c1<<16)|(c2<<8)|c3
cur = {pack(m,0,0,0): 1}
for step in range(t):
nxt = {}
for key, ways in cur.items():
c0=(key>>24)&0xFF; c1=(key>>16)&0xFF; c2=(key>>8)&0xFF; c3=key&0xFF
for x0,x1,x2,x3 in picks:
if x0>c0 or x1>c1 or x2>c2 or x3>c3: continue
nc0=c0-x0; nc1=c1-x1+x0; nc2=c2-x2+x1; nc3=c3-x3+x2
wp = comb[c0][x0]*comb[c1][x1]%MOD*comb[c2][x2]%MOD*comb[c3][x3]%MOD
add = ways * wp % MOD
nk = pack(nc0,nc1,nc2,nc3)
nxt[nk] = (nxt.get(nk, 0) + add) % MOD
cur = nxt
mc = cur.get(pack(0,0,0,0), 0)
aqm = mod_pow(24, m)
ft = 1
for i in range(2, t+1): ft = ft * i % MOD
ift = mod_pow(ft, MOD-2)
return str(mc * aqm % MOD * ift % MOD)
if __name__ == '__main__':
print(solve())
Java
import java.util.HashMap;
import java.util.Map;
public class Euler475 {
private static final long kMod = 1_000_000_007L;
private static long modPow(long base, long exp) {
long out = 1;
base %= kMod;
while (exp > 0) {
if (exp % 2 == 1) {
out = (out * base) % kMod;
}
base = (base * base) % kMod;
exp /= 2;
}
return out;
}
private static int packState(int c0, int c1, int c2, int c3) {
return (c0 << 24) | (c1 << 16) | (c2 << 8) | c3;
}
private static long solve(int musicians) {
int n = musicians / 12;
int m = 3 * n;
int t = 4 * n;
long[][] comb = new long[m + 1][4];
for (int i = 0; i <= m; i++) {
comb[i][0] = 1;
if (i >= 1)
comb[i][1] = i;
if (i >= 2)
comb[i][2] = (long) i * (i - 1) / 2;
if (i >= 3)
comb[i][3] = (long) i * (i - 1) * (i - 2) / 6;
}
int[][] picks = new int[20][4];
int pickCount = 0;
for (int x0 = 0; x0 <= 3; ++x0) {
for (int x1 = 0; x1 + x0 <= 3; ++x1) {
for (int x2 = 0; x2 + x1 + x0 <= 3; ++x2) {
int x3 = 3 - x0 - x1 - x2;
picks[pickCount][0] = x0;
picks[pickCount][1] = x1;
picks[pickCount][2] = x2;
picks[pickCount][3] = x3;
pickCount++;
}
}
}
Map<Integer, Long> cur = new HashMap<>();
Map<Integer, Long> nxt = new HashMap<>();
cur.put(packState(m, 0, 0, 0), 1L);
for (int step = 0; step < t; step++) {
nxt.clear();
for (Map.Entry<Integer, Long> entry : cur.entrySet()) {
int key = entry.getKey();
long waysSoFar = entry.getValue();
int c0 = (key >> 24) & 0xFF;
int c1 = (key >> 16) & 0xFF;
int c2 = (key >> 8) & 0xFF;
int c3 = key & 0xFF;
for (int i = 0; i < pickCount; i++) {
int[] pick = picks[i];
int x0 = pick[0], x1 = pick[1], x2 = pick[2], x3 = pick[3];
if (x0 > c0 || x1 > c1 || x2 > c2 || x3 > c3)
continue;
int nc0 = c0 - x0;
int nc1 = c1 - x1 + x0;
int nc2 = c2 - x2 + x1;
int nc3 = c3 - x3 + x2;
long waysPick = comb[c0][x0];
waysPick = (waysPick * comb[c1][x1]) % kMod;
waysPick = (waysPick * comb[c2][x2]) % kMod;
waysPick = (waysPick * comb[c3][x3]) % kMod;
long add = (waysSoFar * waysPick) % kMod;
int nkey = packState(nc0, nc1, nc2, nc3);
long currentVal = nxt.getOrDefault(nkey, 0L);
nxt.put(nkey, (currentVal + add) % kMod);
}
}
Map<Integer, Long> temp = cur;
cur = nxt;
nxt = temp;
}
long matrixCount = cur.getOrDefault(packState(0, 0, 0, 0), 0L);
long assignQuartetMembers = modPow(24, m);
long factT = 1;
for (int i = 2; i <= t; ++i) {
factT = (factT * i) % kMod;
}
long invFactT = modPow(factT, kMod - 2);
long out = (matrixCount * assignQuartetMembers) % kMod;
out = (out * invFactT) % kMod;
return out;
}
public static void main(String[] args) {
System.out.println(solve(600));
}
}