Problem 696: Mahjong
View on Project EulerProject Euler Problem 696 Solution
EulerSolve provides an optimized solution for Project Euler Problem 696, Mahjong, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary We are given \(s\) suits, each suit has ranks \(1,2,\dots,n\), and each tile type appears at most four times. A valid Mahjong-style hand consists of exactly \(t\) melds together with one pair, where a meld is either a triple of identical tiles or a run of three consecutive ranks inside one suit. The goal is to count all valid hands modulo \(10^9+7\). The implementations evaluate the extremely large case \(W(10^8,10^8,30)\), so any approach based on direct enumeration is hopeless. The entire solution therefore compresses the local suit constraints into a finite automaton, builds one-suit counting polynomials, and only then combines the suits. Mathematical Approach Let \(W(n,s,t)\) denote the number of valid hands with \(n\) ranks per suit, \(s\) suits, and exactly \(t\) melds. The key observation is that suits interact only through the total meld count and the fact that exactly one suit contains the unique pair. Step 1: Split a Suit into Independent Positive Blocks Fix one suit and write its tile multiplicities as \(c_1,\dots,c_n\), where every \(c_i\) lies in \(\{0,1,2,3,4\}\). A zero count breaks all possible runs, so no meld can cross a rank with \(c_i=0\). Therefore a suit decomposes into maximal contiguous blocks of positive counts....
Detailed mathematical approach
Problem Summary
We are given \(s\) suits, each suit has ranks \(1,2,\dots,n\), and each tile type appears at most four times. A valid Mahjong-style hand consists of exactly \(t\) melds together with one pair, where a meld is either a triple of identical tiles or a run of three consecutive ranks inside one suit.
The goal is to count all valid hands modulo \(10^9+7\). The implementations evaluate the extremely large case \(W(10^8,10^8,30)\), so any approach based on direct enumeration is hopeless. The entire solution therefore compresses the local suit constraints into a finite automaton, builds one-suit counting polynomials, and only then combines the suits.
Mathematical Approach
Let \(W(n,s,t)\) denote the number of valid hands with \(n\) ranks per suit, \(s\) suits, and exactly \(t\) melds. The key observation is that suits interact only through the total meld count and the fact that exactly one suit contains the unique pair.
Step 1: Split a Suit into Independent Positive Blocks
Fix one suit and write its tile multiplicities as \(c_1,\dots,c_n\), where every \(c_i\) lies in \(\{0,1,2,3,4\}\). A zero count breaks all possible runs, so no meld can cross a rank with \(c_i=0\).
Therefore a suit decomposes into maximal contiguous blocks of positive counts. Each block can be analyzed independently, and the only remaining global question inside that suit is how many zero gaps are inserted between the blocks. This is why the solution first counts contiguous positive blocks and postpones the placement of zeros until later.
Step 2: Encode One Positive Block with a Small Automaton
Consider a positive block with counts \(d_1,\dots,d_L\), where every \(d_i\in\{1,2,3,4\}\). We process the block from left to right. Before rank \(i\), the local state is \((u,v,\pi)\):
\(u\) is the number of runs started at rank \(i-1\) that still need tiles at ranks \(i\) and \(i+1\).
\(v\) is the number of runs started at rank \(i-2\) that need one final tile at rank \(i\).
\(\pi\in\{0,1\}\) records whether the unique pair has already been used inside the current block.
At rank \(i\), at least \(u+v\) tiles are forced, because those tiles must continue or finish the already open runs. If \(d_i\ge u+v\), let
$$r_i=d_i-u-v.$$
From these remaining tiles we may optionally spend \(2\) tiles on the pair, start some new runs at rank \(i\), and use the rest as triples. If \(\delta\in\{0,1\}\) indicates whether the pair is used at rank \(i\), and \(q_i\) is the number of new runs started there, then the transition is valid exactly when
$$0\le q_i\le r_i-2\delta,\qquad r_i-2\delta-q_i\equiv 0 \pmod 3.$$
The next state becomes
$$ (u,v,\pi)\longrightarrow(q_i,u,\max(\pi,\delta)). $$
Because several choices of \((\delta,q_i)\) may be possible for the same input count \(d_i\), the natural local machine is nondeterministic. The implementations determinize it once and then run ordinary dynamic programming on the reachable deterministic states.
Step 3: Count Contiguous Blocks by Length and Tile Total
Let \(\mathcal{B}_0(L,M)\) be the number of positive blocks of length \(L\) and total tile count \(M\) that end with no unfinished runs and no pair. Let \(\mathcal{B}_1(L,M)\) be the analogous count with exactly one pair used.
These values are obtained by pushing the deterministic automaton through all block lengths up to \(3t+2\) and all tile totals up to \(3t+2\). That cutoff is sufficient because a complete hand with \(t\) melds and one pair contains only \(3t+2\) tiles in total, so no single suit or block can contribute more than that.
The residue classes are forced:
$$\mathcal{B}_0(L,M)\neq 0 \Longrightarrow M\equiv 0 \pmod 3,$$
$$\mathcal{B}_1(L,M)\neq 0 \Longrightarrow M\equiv 2 \pmod 3.$$
Indeed, a block with no pair is built entirely from melds and therefore uses a multiple of three tiles, while a block with one pair uses \(2\) more than a multiple of three.
Step 4: Assemble Several Blocks Inside One Suit
A suit can contain several positive blocks separated by zero-count ranks. Let \(\mathcal{D}_{0,k}(L,M)\) be the number of ways to choose exactly \(k\) blocks whose total positive length is \(L\), total tile count is \(M\), and which contain no pair overall. Let \(\mathcal{D}_{1,k}(L,M)\) be the analogous quantity with exactly one pair overall.
The block convolution is then
$$\mathcal{D}_{0,k+1}(L+L',M+M')\;{+}{=}\;\mathcal{D}_{0,k}(L,M)\,\mathcal{B}_0(L',M'),$$
$$\mathcal{D}_{1,k+1}(L+L',M+M')\;{+}{=}\;\mathcal{D}_{1,k}(L,M)\,\mathcal{B}_0(L',M')+\mathcal{D}_{0,k}(L,M)\,\mathcal{B}_1(L',M').$$
Only one block is allowed to contribute the pair, so the second recurrence has exactly two sources. Also, there can be at most \(t+1\) blocks in one suit, because each nonempty block must contribute at least one group, and a full hand has only \(t\) melds plus one pair.
Step 5: Put Those Blocks Back into the \(n\) Ranks
Suppose a one-suit configuration uses \(k\) positive blocks whose lengths sum to \(L\). The remaining \(n-L\) ranks have count zero. Write the zero gaps as
$$z_0+z_1+\cdots+z_k=n-L,$$
where \(z_0\) is the leading zero run, \(z_k\) is the trailing zero run, and the interior gaps \(z_1,\dots,z_{k-1}\) must each satisfy \(z_i\ge 1\) so that the blocks stay separated.
After setting \(z_i'=z_i-1\) for \(1\le i\le k-1\), we get
$$z_0+z_1'+\cdots+z_{k-1}'+z_k=n-L-k+1,$$
so the number of embeddings is
$$\binom{n-L+1}{k}.$$
This placement factor is one of the main simplifications in the algorithm. For example, if \(n=7\), the chosen blocks have total positive length \(L=4\), and there are \(k=2\) blocks, then the number of placements is \(\binom{7-4+1}{2}=\binom{4}{2}=6\).
Step 6: Convert One Suit into Generating Functions
After multiplying each \(\mathcal{D}_{\varepsilon,k}(L,M)\) by the placement factor, we obtain one-suit totals
$$A_M(n)=\sum_{k,L}\mathcal{D}_{0,k}(L,M)\binom{n-L+1}{k},$$
$$B_M(n)=\sum_{k,L}\mathcal{D}_{1,k}(L,M)\binom{n-L+1}{k}.$$
Now define \(T_u(n)=A_{3u}(n)\) and \(P_u(n)=B_{3u+2}(n)\). In words, \(T_u(n)\) counts one-suit configurations contributing exactly \(u\) melds and no pair, while \(P_u(n)\) counts one-suit configurations contributing exactly \(u\) melds together with the pair.
Introduce the truncated generating functions
$$\mathcal{T}_n(x)=\sum_{u\ge 0}T_u(n)x^u,\qquad \mathcal{P}_n(x)=\sum_{u\ge 0}P_u(n)x^u.$$
Exactly one suit carries the pair, so the full answer is
$$\boxed{W(n,s,t)=s\cdot [x^t]\;\mathcal{P}_n(x)\,\mathcal{T}_n(x)^{\,s-1}\pmod{10^9+7}.}$$
The coefficient of \(x^t\) enforces that the total number of melds across all suits is exactly \(t\).
Worked Example: \(W(4,1,1)=20\)
With one suit, four ranks, and exactly one meld, every legal hand has \(3\cdot 1+2=5\) tiles. There are two cases.
Triple plus pair: choose the triple rank in \(4\) ways and the pair rank in one of the remaining \(3\) ways, giving \(4\cdot 3=12\).
Run plus pair: the only runs are \(123\) and \(234\), so there are \(2\) possible runs. For each run, the pair can be placed at any of the \(4\) ranks, giving \(2\cdot 4=8\).
Therefore
$$W(4,1,1)=12+8=20,$$
which agrees with the small sanity case checked by the implementations.
How the Code Works
The C++, Python, and Java implementations all follow the same pipeline. First they build the reachable deterministic automaton corresponding to the local block rules above. Acceptance is split into two categories: blocks that finish with no pair and blocks that finish with exactly one pair.
Next they run a dynamic program over block length and tile total, but only for positive counts \(1,2,3,4\). Zero-count ranks are deliberately excluded at this stage because zeros are handled later by the explicit placement factor \(\binom{n-L+1}{k}\).
After the contiguous block counts are known, the implementations run a second convolution-style dynamic program that combines many blocks inside one suit, once for the no-pair case and once for the one-pair case. This produces the tables indexed by block count, total positive length, and total tile count.
Then, for each admissible pair \((k,L)\), they multiply by \(\binom{n-L+1}{k}\) modulo \(10^9+7\). Because \(k\le t+1\) is small even when \(n\) is enormous, the binomial coefficient is computed as a short falling product times the modular inverse of \(k!\), rather than by precomputing factorials up to \(n\).
Finally the one-suit coefficients are packed into the truncated polynomials \(\mathcal{T}_n(x)\) and \(\mathcal{P}_n(x)\). The no-pair polynomial is raised to the \((s-1)\)-st power by binary exponentiation with truncation at degree \(t\), multiplied by the pair polynomial, and the coefficient of \(x^t\) is extracted and multiplied by \(s\).
Complexity Analysis
Let \(T=3t+2\). The deterministic automaton has constant size with respect to the input parameters, so the local block precomputation runs in \(O(T^2)\) time and uses \(O(T^2)\) working memory, up to that constant automaton factor.
The dominant part is the block-combination dynamic program. There are \(O(T^2)\) possible block types \((L,M)\), each convolution touches \(O(T^2)\) table cells, and this is repeated for \(O(t)\) block counts. Since \(T=O(t)\), the full precomputation is \(O(t^5)\) time with \(O(t^3)\) stored state.
Once the one-suit tables are built, the final evaluation for given \(n\) and \(s\) takes \(O(t^3+t^2\log s)\) time: \(O(t^3)\) to aggregate the one-suit coefficients and \(O(t^2\log s)\) for truncated polynomial exponentiation. This is exactly why the method remains practical even when \(n\) and \(s\) are both \(10^8\).
Footnotes and References
- Problem page: https://projecteuler.net/problem=696
- Deterministic finite automaton: Wikipedia - Deterministic finite automaton
- Dynamic programming: Wikipedia - Dynamic programming
- Stars and bars: Wikipedia - Stars and bars
- Generating function: Wikipedia - Generating function
- Exponentiation by squaring: Wikipedia - Exponentiation by squaring
Problem 696 source code
C++
#include <algorithm>
#include <array>
#include <cstdint>
#include <iostream>
#include <queue>
#include <thread>
#include <unordered_map>
#include <vector>
namespace {
constexpr int MOD = 1'000'000'007;
long long mod_pow(long long a, long long e) {
long long r = 1 % MOD;
a %= MOD;
while (e > 0) {
if (e & 1LL) r = (r * a) % MOD;
a = (a * a) % MOD;
e >>= 1LL;
}
return r;
}
struct Block {
int len;
int tiles;
int count;
};
struct DFA {
std::vector<std::array<int, 5>> trans;
std::vector<int> accept0;
std::vector<int> accept1;
};
DFA build_dfa() {
auto id = [](int a, int b, int p) { return p * 25 + a * 5 + b; };
const int nfa_states = 50;
std::array<std::array<uint64_t, 5>, nfa_states> nfa{};
for (int p = 0; p <= 1; ++p) {
for (int a = 0; a <= 4; ++a) {
for (int b = 0; b <= 4; ++b) {
const int sid = id(a, b, p);
for (int c = 0; c <= 4; ++c) {
uint64_t mask = 0;
if (c >= a + b) {
const int remaining = c - a - b;
for (int use_pair = 0; use_pair <= 1; ++use_pair) {
if (use_pair && p) continue;
if (use_pair && remaining < 2) continue;
const int rem2 = remaining - 2 * use_pair;
for (int q = 0; q <= 4; ++q) {
if (q > rem2) break;
if ((rem2 - q) % 3 != 0) continue;
const int ns = id(q, a, p || use_pair);
mask |= 1ULL << ns;
}
}
}
nfa[sid][c] = mask;
}
}
}
}
const uint64_t start_mask = 1ULL << id(0, 0, 0);
std::unordered_map<uint64_t, int> index;
std::vector<uint64_t> masks;
std::queue<uint64_t> q;
index[start_mask] = 0;
masks.push_back(start_mask);
q.push(start_mask);
std::vector<std::array<int, 5>> trans;
while (!q.empty()) {
const uint64_t mask = q.front();
q.pop();
const int cur = index[mask];
if (cur >= static_cast<int>(trans.size())) {
trans.push_back({0, 0, 0, 0, 0});
}
for (int c = 0; c <= 4; ++c) {
uint64_t next_mask = 0;
uint64_t tmp = mask;
while (tmp) {
const int i = __builtin_ctzll(tmp);
tmp &= tmp - 1;
next_mask |= nfa[i][c];
}
auto it = index.find(next_mask);
int nxt;
if (it == index.end()) {
nxt = static_cast<int>(index.size());
index[next_mask] = nxt;
masks.push_back(next_mask);
q.push(next_mask);
} else {
nxt = it->second;
}
trans[cur][c] = nxt;
}
}
const uint64_t acc0_bit = 1ULL << id(0, 0, 0);
const uint64_t acc1_bit = 1ULL << id(0, 0, 1);
std::vector<int> accept0;
std::vector<int> accept1;
for (int i = 0; i < static_cast<int>(masks.size()); ++i) {
if (masks[i] & acc0_bit) accept0.push_back(i);
if (masks[i] & acc1_bit) accept1.push_back(i);
}
DFA dfa;
dfa.trans = std::move(trans);
dfa.accept0 = std::move(accept0);
dfa.accept1 = std::move(accept1);
return dfa;
}
class MahjongCounter {
public:
MahjongCounter(int max_t, unsigned threads)
: max_t_(max_t),
max_tiles_(3 * max_t + 2),
max_len_(3 * max_t + 2),
max_blocks_(max_t + 1),
stride_(max_tiles_ + 1),
threads_(std::max(1u, threads)),
dfa_(build_dfa()) {
build_factorials();
build_block_counts();
build_block_dp();
}
long long solve(long long n, long long s, int t) const {
if (s == 0) return 0;
if (t > max_t_) return 0;
const int tiles_max = 3 * t + 2;
const int len_max = tiles_max;
const int blocks_max = t + 1;
std::vector<int> A(tiles_max + 1, 0);
std::vector<int> B(tiles_max + 1, 0);
for (int k = 0; k <= blocks_max; ++k) {
for (int L = 0; L <= len_max; ++L) {
const long long N = n - L + 1;
if (N < k) continue;
const int comb = comb_small(N, k);
if (comb == 0) continue;
const int base = L * stride_;
const std::vector<int>& dp0 = dp0_[k];
const std::vector<int>& dp1 = dp1_[k];
for (int M = 0; M <= tiles_max; ++M) {
const int v0 = dp0[base + M];
if (v0) {
A[M] = (A[M] + static_cast<long long>(v0) * comb) % MOD;
}
const int v1 = dp1[base + M];
if (v1) {
B[M] = (B[M] + static_cast<long long>(v1) * comb) % MOD;
}
}
}
}
std::vector<int> G(t + 1, 0);
std::vector<int> H(t + 1, 0);
for (int k = 0; k <= t; ++k) {
const int idxG = 3 * k;
const int idxH = 3 * k + 2;
if (idxG <= tiles_max) G[k] = A[idxG];
if (idxH <= tiles_max) H[k] = B[idxH];
}
std::vector<int> Gpow = poly_pow(G, s - 1, t);
long long total = 0;
for (int i = 0; i <= t; ++i) {
total = (total + static_cast<long long>(H[i]) * Gpow[t - i]) % MOD;
}
total = total * (s % MOD) % MOD;
return total;
}
private:
int max_t_;
int max_tiles_;
int max_len_;
int max_blocks_;
int stride_;
unsigned threads_;
DFA dfa_;
std::vector<int> fact_;
std::vector<int> inv_fact_;
std::vector<std::vector<int>> block0_;
std::vector<std::vector<int>> block1_;
std::vector<Block> triple_blocks_;
std::vector<Block> pair_blocks_;
std::vector<std::vector<int>> dp0_;
std::vector<std::vector<int>> dp1_;
void build_factorials() {
fact_.assign(max_blocks_ + 1, 1);
inv_fact_.assign(max_blocks_ + 1, 1);
for (int i = 1; i <= max_blocks_; ++i) {
fact_[i] = static_cast<long long>(fact_[i - 1]) * i % MOD;
}
inv_fact_[max_blocks_] = static_cast<int>(mod_pow(fact_[max_blocks_], MOD - 2));
for (int i = max_blocks_; i >= 1; --i) {
inv_fact_[i - 1] = static_cast<long long>(inv_fact_[i]) * i % MOD;
}
}
int comb_small(long long N, int k) const {
if (k < 0 || N < k) return 0;
long long res = 1;
for (int i = 0; i < k; ++i) {
long long term = (N - i) % MOD;
if (term < 0) term += MOD;
res = res * term % MOD;
}
res = res * inv_fact_[k] % MOD;
return static_cast<int>(res);
}
void build_block_counts() {
const int S = static_cast<int>(dfa_.trans.size());
const int size = S * stride_;
std::vector<int> cur(size, 0);
std::vector<int> next(size, 0);
block0_.assign(max_len_ + 1, std::vector<int>(max_tiles_ + 1, 0));
block1_.assign(max_len_ + 1, std::vector<int>(max_tiles_ + 1, 0));
cur[0] = 1;
for (int L = 1; L <= max_len_; ++L) {
std::fill(next.begin(), next.end(), 0);
for (int s = 0; s < S; ++s) {
const int base = s * stride_;
for (int m = 0; m <= max_tiles_; ++m) {
const int val = cur[base + m];
if (!val) continue;
for (int c = 1; c <= 4; ++c) {
const int nm = m + c;
if (nm > max_tiles_) break;
const int ns = dfa_.trans[s][c];
const int idx = ns * stride_ + nm;
int v = next[idx] + val;
if (v >= MOD) v -= MOD;
next[idx] = v;
}
}
}
std::vector<int> sum0(max_tiles_ + 1, 0);
std::vector<int> sum1(max_tiles_ + 1, 0);
for (int s : dfa_.accept0) {
const int base = s * stride_;
for (int m = 0; m <= max_tiles_; ++m) {
int v = sum0[m] + next[base + m];
if (v >= MOD) v -= MOD;
sum0[m] = v;
}
}
for (int s : dfa_.accept1) {
const int base = s * stride_;
for (int m = 0; m <= max_tiles_; ++m) {
int v = sum1[m] + next[base + m];
if (v >= MOD) v -= MOD;
sum1[m] = v;
}
}
block0_[L] = std::move(sum0);
block1_[L] = std::move(sum1);
cur.swap(next);
}
for (int L = 1; L <= max_len_; ++L) {
for (int M = 0; M <= max_tiles_; ++M) {
const int c0 = block0_[L][M];
if (c0 && (M % 3 == 0)) {
triple_blocks_.push_back({L, M, c0});
}
const int c1 = block1_[L][M];
if (c1 && (M % 3 == 2)) {
pair_blocks_.push_back({L, M, c1});
}
}
}
}
void convolve_add(const std::vector<int>& dp,
const std::vector<Block>& blocks,
std::vector<int>& out) const {
if (blocks.empty()) return;
if (threads_ <= 1 || blocks.size() < 200) {
for (const auto& b : blocks) {
const int Lb = b.len;
const int Mb = b.tiles;
const int cnt = b.count;
const int Llimit = max_len_ - Lb;
const int Mlimit = max_tiles_ - Mb;
for (int L = 0; L <= Llimit; ++L) {
const int base = L * stride_;
const int out_base = (L + Lb) * stride_ + Mb;
for (int M = 0; M <= Mlimit; ++M) {
const int val = dp[base + M];
if (!val) continue;
const int idx = out_base + M;
out[idx] = (out[idx] + static_cast<long long>(val) * cnt) % MOD;
}
}
}
return;
}
const size_t total = out.size();
const unsigned tcount = std::min<unsigned>(threads_, blocks.size());
std::vector<std::vector<int>> locals(tcount, std::vector<int>(total, 0));
std::vector<std::thread> workers;
workers.reserve(tcount);
for (unsigned t = 0; t < tcount; ++t) {
const size_t start = t * blocks.size() / tcount;
const size_t end = (t + 1) * blocks.size() / tcount;
workers.emplace_back([&, start, end, t]() {
auto& local = locals[t];
for (size_t i = start; i < end; ++i) {
const auto& b = blocks[i];
const int Lb = b.len;
const int Mb = b.tiles;
const int cnt = b.count;
const int Llimit = max_len_ - Lb;
const int Mlimit = max_tiles_ - Mb;
for (int L = 0; L <= Llimit; ++L) {
const int base = L * stride_;
const int out_base = (L + Lb) * stride_ + Mb;
for (int M = 0; M <= Mlimit; ++M) {
const int val = dp[base + M];
if (!val) continue;
const int idx = out_base + M;
local[idx] = (local[idx] + static_cast<long long>(val) * cnt) % MOD;
}
}
}
});
}
for (auto& th : workers) th.join();
for (size_t i = 0; i < total; ++i) {
long long sum = out[i];
for (unsigned t = 0; t < tcount; ++t) {
sum += locals[t][i];
}
out[i] = static_cast<int>(sum % MOD);
}
}
void build_block_dp() {
const int size = (max_len_ + 1) * stride_;
dp0_.assign(max_blocks_ + 1, std::vector<int>(size, 0));
dp1_.assign(max_blocks_ + 1, std::vector<int>(size, 0));
dp0_[0][0] = 1;
for (int k = 0; k < max_blocks_; ++k) {
std::vector<int> next0(size, 0);
std::vector<int> next1(size, 0);
convolve_add(dp0_[k], triple_blocks_, next0);
convolve_add(dp1_[k], triple_blocks_, next1);
convolve_add(dp0_[k], pair_blocks_, next1);
dp0_[k + 1] = std::move(next0);
dp1_[k + 1] = std::move(next1);
}
}
static std::vector<int> poly_mul(const std::vector<int>& a,
const std::vector<int>& b,
int t) {
std::vector<int> res(t + 1, 0);
const int n = std::min<int>(t, static_cast<int>(a.size()) - 1);
const int m = std::min<int>(t, static_cast<int>(b.size()) - 1);
for (int i = 0; i <= n; ++i) {
if (!a[i]) continue;
for (int j = 0; j <= m && i + j <= t; ++j) {
if (!b[j]) continue;
res[i + j] = (res[i + j] + static_cast<long long>(a[i]) * b[j]) % MOD;
}
}
return res;
}
static std::vector<int> poly_pow(std::vector<int> base, long long exp, int t) {
std::vector<int> res(t + 1, 0);
res[0] = 1;
while (exp > 0) {
if (exp & 1LL) res = poly_mul(res, base, t);
exp >>= 1LL;
if (exp) base = poly_mul(base, base, t);
}
return res;
}
};
} // namespace
int main() {
const int max_t = 30;
unsigned threads = std::thread::hardware_concurrency();
if (threads == 0) threads = 1;
threads = std::min(threads, 8u);
MahjongCounter counter(max_t, threads);
bool ok = true;
auto check = [&](long long n, long long s, int t, long long expected) {
const long long got = counter.solve(n, s, t);
if (got != expected) {
std::cerr << "[CHECK FAILED] w(" << n << ", " << s << ", " << t << ") = "
<< got << ", expected " << expected << "\n";
ok = false;
}
};
check(4, 1, 1, 20);
check(9, 1, 4, 13259);
check(9, 3, 4, 5237550);
check(1000, 1000, 5, 107662178);
if (!ok) return 1;
const long long n = 100000000LL;
const long long s = 100000000LL;
const int t = 30;
const long long ans = counter.solve(n, s, t);
std::cout << ans % MOD << "\n";
return 0;
}
Python
import sys
import collections
MOD = 1000000007
def mod_pow(a, e):
r = 1
a %= MOD
while e > 0:
if e & 1: r = (r * a) % MOD
a = (a * a) % MOD
e >>= 1
return r
def build_dfa():
def get_id(a, b, p):
return p * 25 + a * 5 + b
nfa_states = 50
nfa = [[0] * 5 for _ in range(nfa_states)]
for p in range(2):
for a in range(5):
for b in range(5):
sid = get_id(a, b, p)
for c in range(5):
mask = 0
if c >= a + b:
remaining = c - a - b
for use_pair in range(2):
if use_pair and p: continue
if use_pair and remaining < 2: continue
rem2 = remaining - 2 * use_pair
for q in range(5):
if q > rem2: break
if (rem2 - q) % 3 != 0: continue
ns = get_id(q, a, p or use_pair)
mask |= (1 << ns)
nfa[sid][c] = mask
start_mask = 1 << get_id(0, 0, 0)
index_map = {start_mask: 0}
masks = [start_mask]
q = collections.deque([start_mask])
trans = []
while q:
mask = q.popleft()
cur = index_map[mask]
if cur >= len(trans):
trans.append([0] * 5)
for c in range(5):
next_mask = 0
tmp = mask
while tmp:
i = (tmp & -tmp).bit_length() - 1
tmp &= tmp - 1
next_mask |= nfa[i][c]
if next_mask not in index_map:
nxt = len(index_map)
index_map[next_mask] = nxt
masks.append(next_mask)
q.append(next_mask)
else:
nxt = index_map[next_mask]
trans[cur][c] = nxt
acc0_bit = 1 << get_id(0, 0, 0)
acc1_bit = 1 << get_id(0, 0, 1)
accept0 = [i for i, m in enumerate(masks) if (m & acc0_bit)]
accept1 = [i for i, m in enumerate(masks) if (m & acc1_bit)]
return trans, accept0, accept1
class MahjongCounter:
def __init__(self, max_t):
self.max_t = max_t
self.max_tiles = 3 * max_t + 2
self.max_len = 3 * max_t + 2
self.max_blocks = max_t + 1
self.stride = self.max_tiles + 1
self.trans, self.accept0, self.accept1 = build_dfa()
self.fact = [1] * (self.max_blocks + 1)
self.inv_fact = [1] * (self.max_blocks + 1)
for i in range(1, self.max_blocks + 1):
self.fact[i] = (self.fact[i - 1] * i) % MOD
self.inv_fact[self.max_blocks] = mod_pow(self.fact[self.max_blocks], MOD - 2)
for i in range(self.max_blocks, 0, -1):
self.inv_fact[i - 1] = (self.inv_fact[i] * i) % MOD
self.build_block_counts()
self.build_block_dp()
def comb_small(self, N, k):
if k < 0 or N < k: return 0
res = 1
for i in range(k):
term = (N - i) % MOD
if term < 0: term += MOD
res = (res * term) % MOD
res = (res * self.inv_fact[k]) % MOD
return res
def build_block_counts(self):
S = len(self.trans)
size = S * self.stride
cur = [0] * size
cur[0] = 1
self.triple_blocks = []
self.pair_blocks = []
for L in range(1, self.max_len + 1):
nxt = [0] * size
for s in range(S):
base = s * self.stride
for m in range(self.max_tiles + 1):
val = cur[base + m]
if not val: continue
for c in range(1, 5):
nm = m + c
if nm > self.max_tiles: break
ns = self.trans[s][c]
idx = ns * self.stride + nm
nxt[idx] = (nxt[idx] + val) % MOD
sum0 = [0] * (self.max_tiles + 1)
sum1 = [0] * (self.max_tiles + 1)
for s in self.accept0:
base = s * self.stride
for m in range(self.max_tiles + 1):
sum0[m] = (sum0[m] + nxt[base + m]) % MOD
for s in self.accept1:
base = s * self.stride
for m in range(self.max_tiles + 1):
sum1[m] = (sum1[m] + nxt[base + m]) % MOD
for M in range(self.max_tiles + 1):
if sum0[M] and M % 3 == 0:
self.triple_blocks.append((L, M, sum0[M]))
if sum1[M] and M % 3 == 2:
self.pair_blocks.append((L, M, sum1[M]))
cur = nxt
def convolve_add(self, dp, blocks, out):
if not blocks: return
for Lb, Mb, cnt in blocks:
Llimit = self.max_len - Lb
Mlimit = self.max_tiles - Mb
for L in range(Llimit + 1):
base = L * self.stride
out_base = (L + Lb) * self.stride + Mb
for M in range(Mlimit + 1):
val = dp[base + M]
if not val: continue
idx = out_base + M
out[idx] = (out[idx] + val * cnt) % MOD
def build_block_dp(self):
size = (self.max_len + 1) * self.stride
self.dp0 = [[0] * size for _ in range(self.max_blocks + 1)]
self.dp1 = [[0] * size for _ in range(self.max_blocks + 1)]
self.dp0[0][0] = 1
for k in range(self.max_blocks):
self.convolve_add(self.dp0[k], self.triple_blocks, self.dp0[k + 1])
self.convolve_add(self.dp1[k], self.triple_blocks, self.dp1[k + 1])
self.convolve_add(self.dp0[k], self.pair_blocks, self.dp1[k + 1])
def poly_mul(self, a, b, t):
res = [0] * (t + 1)
for i in range(min(t + 1, len(a))):
if not a[i]: continue
for j in range(min(t + 1 - i, len(b))):
if not b[j]: continue
res[i + j] = (res[i + j] + a[i] * b[j]) % MOD
return res
def poly_pow(self, base, exp, t):
res = [0] * (t + 1)
res[0] = 1
while exp > 0:
if exp & 1: res = self.poly_mul(res, base, t)
exp >>= 1
if exp: base = self.poly_mul(base, base, t)
return res
def solve(self, n, s, t):
if s == 0 or t > self.max_t: return 0
A = [0] * (self.max_tiles + 1)
B = [0] * (self.max_tiles + 1)
for k in range(self.max_blocks + 1):
for L in range(self.max_len + 1):
N_val = n - L + 1
if N_val < k: continue
comb = self.comb_small(N_val, k)
if comb == 0: continue
base = L * self.stride
dp0_k = self.dp0[k]
dp1_k = self.dp1[k]
for M in range(self.max_tiles + 1):
if dp0_k[base + M]:
A[M] = (A[M] + dp0_k[base + M] * comb) % MOD
if dp1_k[base + M]:
B[M] = (B[M] + dp1_k[base + M] * comb) % MOD
G = [0] * (t + 1)
H = [0] * (t + 1)
for k in range(t + 1):
idxG = 3 * k
idxH = 3 * k + 2
if idxG <= self.max_tiles: G[k] = A[idxG]
if idxH <= self.max_tiles: H[k] = B[idxH]
Gpow = self.poly_pow(G, s - 1, t)
total = 0
for i in range(t + 1):
total = (total + H[i] * Gpow[t - i]) % MOD
total = (total * (s % MOD)) % MOD
return total
def solve():
n = 100000000
s = 100000000
t = 30
counter = MahjongCounter(t)
return str(counter.solve(n, s, t))
if __name__ == '__main__':
print(solve())
Java
import java.util.*;
public class Euler696 {
static final long MOD = 1000000007L;
static long modPow(long a, long e) {
long r = 1;
a %= MOD;
while (e > 0) {
if ((e & 1) != 0)
r = (r * a) % MOD;
a = (a * a) % MOD;
e >>= 1;
}
return r;
}
static class Block {
int len, tiles, count;
Block(int l, int t, int c) {
len = l;
tiles = t;
count = c;
}
}
static class DFA {
int[][] trans;
int[] accept0;
int[] accept1;
}
static DFA buildDFA() {
int nfaStates = 50;
long[][] nfa = new long[nfaStates][5];
for (int p = 0; p <= 1; ++p) {
for (int a = 0; a <= 4; ++a) {
for (int b = 0; b <= 4; ++b) {
int sid = p * 25 + a * 5 + b;
for (int c = 0; c <= 4; ++c) {
long mask = 0;
if (c >= a + b) {
int remaining = c - a - b;
for (int usePair = 0; usePair <= 1; ++usePair) {
if (usePair == 1 && p == 1)
continue;
if (usePair == 1 && remaining < 2)
continue;
int rem2 = remaining - 2 * usePair;
for (int q = 0; q <= 4; ++q) {
if (q > rem2)
break;
if ((rem2 - q) % 3 != 0)
continue;
int ns = (p | usePair) * 25 + q * 5 + a;
mask |= (1L << ns);
}
}
}
nfa[sid][c] = mask;
}
}
}
}
long startMask = 1L << 0;
HashMap<Long, Integer> index = new HashMap<>();
List<Long> masks = new ArrayList<>();
Queue<Long> q = new LinkedList<>();
index.put(startMask, 0);
masks.add(startMask);
q.add(startMask);
List<int[]> trans = new ArrayList<>();
while (!q.isEmpty()) {
long mask = q.poll();
int cur = index.get(mask);
if (cur >= trans.size()) {
trans.add(new int[5]);
}
for (int c = 0; c <= 4; ++c) {
long nextMask = 0;
long tmp = mask;
while (tmp != 0) {
int i = Long.numberOfTrailingZeros(tmp);
tmp &= tmp - 1;
nextMask |= nfa[i][c];
}
int nxt;
if (!index.containsKey(nextMask)) {
nxt = index.size();
index.put(nextMask, nxt);
masks.add(nextMask);
q.add(nextMask);
} else {
nxt = index.get(nextMask);
}
trans.get(cur)[c] = nxt;
}
}
long acc0Bit = 1L << 0;
long acc1Bit = 1L << 25;
List<Integer> accept0List = new ArrayList<>();
List<Integer> accept1List = new ArrayList<>();
for (int i = 0; i < masks.size(); ++i) {
if ((masks.get(i) & acc0Bit) != 0)
accept0List.add(i);
if ((masks.get(i) & acc1Bit) != 0)
accept1List.add(i);
}
DFA dfa = new DFA();
dfa.trans = trans.toArray(new int[0][]);
dfa.accept0 = accept0List.stream().mapToInt(i -> i).toArray();
dfa.accept1 = accept1List.stream().mapToInt(i -> i).toArray();
return dfa;
}
static class MahjongCounter {
int maxT, maxTiles, maxLen, maxBlocks, stride;
DFA dfa;
int[] fact, invFact;
List<Block> tripleBlocks = new ArrayList<>();
List<Block> pairBlocks = new ArrayList<>();
int[][] dp0;
int[][] dp1;
MahjongCounter(int t) {
maxT = t;
maxTiles = 3 * t + 2;
maxLen = 3 * t + 2;
maxBlocks = t + 1;
stride = maxTiles + 1;
dfa = buildDFA();
buildFactorials();
buildBlockCounts();
buildBlockDP();
}
void buildFactorials() {
fact = new int[maxBlocks + 1];
invFact = new int[maxBlocks + 1];
fact[0] = 1;
for (int i = 1; i <= maxBlocks; ++i) {
fact[i] = (int) ((fact[i - 1] * (long) i) % MOD);
}
invFact[maxBlocks] = (int) modPow(fact[maxBlocks], MOD - 2);
for (int i = maxBlocks; i >= 1; --i) {
invFact[i - 1] = (int) ((invFact[i] * (long) i) % MOD);
}
}
int combSmall(long N, int k) {
if (k < 0 || N < k)
return 0;
long res = 1;
for (int i = 0; i < k; ++i) {
long term = (N - i) % MOD;
if (term < 0)
term += MOD;
res = (res * term) % MOD;
}
res = (res * invFact[k]) % MOD;
return (int) res;
}
void buildBlockCounts() {
int S = dfa.trans.length;
int size = S * stride;
int[] cur = new int[size];
int[] next = new int[size];
cur[0] = 1;
for (int L = 1; L <= maxLen; ++L) {
Arrays.fill(next, 0);
for (int s = 0; s < S; ++s) {
int base = s * stride;
for (int m = 0; m <= maxTiles; ++m) {
int val = cur[base + m];
if (val == 0)
continue;
for (int c = 1; c <= 4; ++c) {
int nm = m + c;
if (nm > maxTiles)
break;
int ns = dfa.trans[s][c];
int idx = ns * stride + nm;
int v = next[idx] + val;
if (v >= MOD)
v -= MOD;
next[idx] = v;
}
}
}
int[] sum0 = new int[maxTiles + 1];
int[] sum1 = new int[maxTiles + 1];
for (int s : dfa.accept0) {
int base = s * stride;
for (int m = 0; m <= maxTiles; ++m) {
int v = sum0[m] + next[base + m];
if (v >= MOD)
v -= MOD;
sum0[m] = v;
}
}
for (int s : dfa.accept1) {
int base = s * stride;
for (int m = 0; m <= maxTiles; ++m) {
int v = sum1[m] + next[base + m];
if (v >= MOD)
v -= MOD;
sum1[m] = v;
}
}
for (int M = 0; M <= maxTiles; ++M) {
if (sum0[M] > 0 && M % 3 == 0)
tripleBlocks.add(new Block(L, M, sum0[M]));
if (sum1[M] > 0 && M % 3 == 2)
pairBlocks.add(new Block(L, M, sum1[M]));
}
int[] temp = cur;
cur = next;
next = temp;
}
}
void convolveAdd(int[] dp, List<Block> blocks, int[] out) {
for (Block b : blocks) {
int Llimit = maxLen - b.len;
int Mlimit = maxTiles - b.tiles;
for (int L = 0; L <= Llimit; ++L) {
int base = L * stride;
int outBase = (L + b.len) * stride + b.tiles;
for (int M = 0; M <= Mlimit; ++M) {
int val = dp[base + M];
if (val == 0)
continue;
int idx = outBase + M;
out[idx] = (int) ((out[idx] + (long) val * b.count) % MOD);
}
}
}
}
void buildBlockDP() {
int size = (maxLen + 1) * stride;
dp0 = new int[maxBlocks + 1][size];
dp1 = new int[maxBlocks + 1][size];
dp0[0][0] = 1;
for (int k = 0; k < maxBlocks; ++k) {
convolveAdd(dp0[k], tripleBlocks, dp0[k + 1]);
convolveAdd(dp1[k], tripleBlocks, dp1[k + 1]);
convolveAdd(dp0[k], pairBlocks, dp1[k + 1]);
}
}
long[] polyMul(long[] a, long[] b, int t) {
long[] res = new long[t + 1];
for (int i = 0; i <= Math.min(t, a.length - 1); ++i) {
if (a[i] == 0)
continue;
for (int j = 0; j <= Math.min(t - i, b.length - 1); ++j) {
if (b[j] == 0)
continue;
res[i + j] = (res[i + j] + a[i] * b[j]) % MOD;
}
}
return res;
}
long[] polyPow(long[] base, long exp, int t) {
long[] res = new long[t + 1];
res[0] = 1;
while (exp > 0) {
if ((exp & 1) != 0)
res = polyMul(res, base, t);
exp >>= 1;
if (exp > 0)
base = polyMul(base, base, t);
}
return res;
}
long solve(long n, long s, int t) {
if (s == 0 || t > maxT)
return 0;
long[] A = new long[maxTiles + 1];
long[] B = new long[maxTiles + 1];
for (int k = 0; k <= maxBlocks; ++k) {
for (int L = 0; L <= maxLen; ++L) {
long N = n - L + 1;
if (N < k)
continue;
int comb = combSmall(N, k);
if (comb == 0)
continue;
int base = L * stride;
int[] dp0K = dp0[k];
int[] dp1K = dp1[k];
for (int M = 0; M <= maxTiles; ++M) {
int v0 = dp0K[base + M];
if (v0 > 0)
A[M] = (A[M] + (long) v0 * comb) % MOD;
int v1 = dp1K[base + M];
if (v1 > 0)
B[M] = (B[M] + (long) v1 * comb) % MOD;
}
}
}
long[] G = new long[t + 1];
long[] H = new long[t + 1];
for (int k = 0; k <= t; ++k) {
int idxG = 3 * k;
int idxH = 3 * k + 2;
if (idxG <= maxTiles)
G[k] = A[idxG];
if (idxH <= maxTiles)
H[k] = B[idxH];
}
long[] Gpow = polyPow(G, s - 1, t);
long total = 0;
for (int i = 0; i <= t; ++i) {
total = (total + H[i] * Gpow[t - i]) % MOD;
}
total = (total * (s % MOD)) % MOD;
return total;
}
}
public static String solve() {
long n = 100000000L;
long s = 100000000L;
int t = 30;
MahjongCounter mc = new MahjongCounter(t);
return Long.toString(mc.solve(n, s, t));
}
public static void main(String[] args) {
System.out.println(solve());
}
}