Problem 677: Coloured Graphs
View on Project EulerProject Euler Problem 677 Solution
EulerSolve provides an optimized solution for Project Euler Problem 677, Coloured Graphs, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary The solution counts admissible unlabelled coloured trees on \(n\) vertices and asks for \(g(10000)\) modulo $$10^9+7.$$ The colour rules encoded by the implementations are: red vertices may have degree at most \(4\), blue vertices degree at most \(3\), yellow vertices degree at most \(3\), and an edge joining two yellow vertices is forbidden. The first checkpoints are $$g(2)=5,\qquad g(3)=15,\qquad g(4)=57,\qquad g(10)=710249.$$ A direct isomorphism test over all coloured trees of size \(10000\) is hopeless, so the program instead uses ordinary generating functions, cycle-index formulas for multisets, and the dissymmetry theorem for trees. Mathematical Approach The counting is done in three stages: first count planted trees, then count vertex-rooted trees, and finally recover the unrooted objects. Step 1: Planted Trees and Colour Constraints Let \(R(x)\), \(B(x)\), and \(Y(x)\) be the ordinary generating functions for planted trees whose distinguished root vertex is red, blue, or yellow. A planted tree has one distinguished parent edge above the root, so one incidence of the root is already used. Define the aggregate series $$A(x)=R(x)+B(x)+Y(x),\qquad N(x)=R(x)+B(x).$$ Here \(A(x)\) is the class of all planted trees, while \(N(x)\) is the class of planted trees whose root is not yellow....
Detailed mathematical approach
Problem Summary
The solution counts admissible unlabelled coloured trees on \(n\) vertices and asks for \(g(10000)\) modulo
$$10^9+7.$$
The colour rules encoded by the implementations are:
red vertices may have degree at most \(4\), blue vertices degree at most \(3\), yellow vertices degree at most \(3\), and an edge joining two yellow vertices is forbidden. The first checkpoints are
$$g(2)=5,\qquad g(3)=15,\qquad g(4)=57,\qquad g(10)=710249.$$
A direct isomorphism test over all coloured trees of size \(10000\) is hopeless, so the program instead uses ordinary generating functions, cycle-index formulas for multisets, and the dissymmetry theorem for trees.
Mathematical Approach
The counting is done in three stages: first count planted trees, then count vertex-rooted trees, and finally recover the unrooted objects.
Step 1: Planted Trees and Colour Constraints
Let \(R(x)\), \(B(x)\), and \(Y(x)\) be the ordinary generating functions for planted trees whose distinguished root vertex is red, blue, or yellow. A planted tree has one distinguished parent edge above the root, so one incidence of the root is already used.
Define the aggregate series
$$A(x)=R(x)+B(x)+Y(x),\qquad N(x)=R(x)+B(x).$$
Here \(A(x)\) is the class of all planted trees, while \(N(x)\) is the class of planted trees whose root is not yellow. The second series is needed because a yellow root is not allowed to touch a yellow neighbour.
Since the parent edge already occupies one slot, the root of a planted tree can accept only
red: \(0,1,2,3\) children, blue: \(0,1,2\) children, yellow: \(0,1,2\) non-yellow children.
Step 2: Multisets via Cycle Indices
When the root receives an unordered collection of child subtrees, we must count multisets rather than ordered tuples. For ordinary generating functions, the cycle-index formulas for the symmetric groups give
$$\Phi_1(F)=F(x),$$
$$\Phi_2(F)=\frac{F(x)^2+F(x^2)}{2},$$
$$\Phi_3(F)=\frac{F(x)^3+3F(x)F(x^2)+2F(x^3)}{6},$$
$$\Phi_4(F)=\frac{F(x)^4+6F(x)^2F(x^2)+3F(x^2)^2+8F(x)F(x^3)+6F(x^4)}{24}.$$
These are exactly the multiset operators that appear in the program.
Therefore the planted-tree equations are
$$R(x)=x\left(1+\Phi_1(A)+\Phi_2(A)+\Phi_3(A)\right),$$
$$B(x)=x\left(1+\Phi_1(A)+\Phi_2(A)\right),$$
$$Y(x)=x\left(1+\Phi_1(N)+\Phi_2(N)\right).$$
The factor \(x\) accounts for the root itself, and the constant \(1\) means that the root may have no children.
Step 3: Vertex-Rooted Trees
To count trees rooted at a vertex rather than at a parent edge, the root regains its full degree budget. Let \(V_R(x)\), \(V_B(x)\), and \(V_Y(x)\) denote the generating functions for vertex-rooted admissible trees. Then
$$V_R(x)=x\left(1+\Phi_1(A)+\Phi_2(A)+\Phi_3(A)+\Phi_4(A)\right),$$
$$V_B(x)=x\left(1+\Phi_1(A)+\Phi_2(A)+\Phi_3(A)\right),$$
$$V_Y(x)=x\left(1+\Phi_1(N)+\Phi_2(N)+\Phi_3(N)\right).$$
Adding the three root colours gives the full vertex-rooted series
$$V(x)=V_R(x)+V_B(x)+V_Y(x).$$
This matches the second stage of the implementations.
Step 4: Edge-Rooted Trees and the Yellow-Yellow Exclusion
If we root a tree at a directed edge, the edge splits the tree into an ordered pair of planted trees. Without the colour restriction this would give \(A(x)^2\). The only forbidden case is when both exposed roots are yellow, so the directed edge-rooted series is
$$E_{\mathrm{dir}}(x)=A(x)^2-Y(x)^2.$$
For an undirected edge root, the two sides are an unordered pair. Burnside's lemma gives
$$E_{\mathrm{und}}(x)=\frac{E_{\mathrm{dir}}(x)+A(x^2)-Y(x^2)}{2}.$$
The extra terms \(A(x^2)\) and \(Y(x^2)\) represent the fixed points of swapping the two sides: both halves must be identical, and the yellow-yellow case is still excluded.
Step 5: Dissymmetry Theorem
For trees, the classical dissymmetry theorem says that
$$\text{unrooted}=\text{vertex-rooted}+\text{edge-rooted}-\text{directed-edge-rooted}.$$
Translated into generating functions,
$$G(x)=V(x)+E_{\mathrm{und}}(x)-E_{\mathrm{dir}}(x).$$
The desired sequence is therefore
$$g(n)=\left[x^n\right]G(x).$$
This identity is exactly what the final correction stage computes.
Worked Example: Why \(g(3)=15\)
For \(n=3\), the only tree shape is a path of length \(2\). Up to isomorphism, the middle vertex determines the adjacency pattern.
If the middle vertex is red, the two leaves form an unordered multiset chosen from \(\{R,B,Y\}\). There are
$$\binom{3+2-1}{2}=6$$
such multisets.
If the middle vertex is blue, the same count applies, again giving \(6\).
If the middle vertex is yellow, each neighbour must be non-yellow, so the two leaves form an unordered multiset from \(\{R,B\}\). That gives
$$\binom{2+2-1}{2}=3.$$
Hence
$$g(3)=6+6+3=15,$$
which agrees with the validation values used by the implementation.
How the Code Works
The C++, Python, and Java implementations all follow the same pipeline. First they build coefficient tables for the planted series up to size \(10000\). When constructing size \(n\), they read multiset contributions at size \(n-1\) because the root contributes one vertex.
During that first pass, the implementation incrementally maintains square convolutions of the aggregated planted series. Those precomputed quadratic terms are then reused to obtain the degree-\(2\) and degree-\(3\) multiset contributions efficiently.
In the second pass, the implementation constructs shifted copies corresponding to \(F(x^2)\), \(F(x^3)\), and \(F(x^4)\), together with the cubic and quartic convolution terms required by the cycle-index formulas. That produces the full vertex-rooted counts.
Finally it forms the directed edge-rooted and undirected edge-rooted series, applies
$$G(x)=V(x)+E_{\mathrm{und}}(x)-E_{\mathrm{dir}}(x),$$
and extracts the coefficient of \(x^{10000}\) modulo \(10^9+7\). The C++ implementation parallelizes some of the largest convolutions, but the mathematical formula is identical in all three languages.
Complexity Analysis
Let \(N=10000\). The dominant work comes from quadratic convolutions and from the incremental updates that build the square series. Both stages are \(O(N^2)\) time overall. Only a fixed number of coefficient arrays of length \(N+1\) are stored, so the memory usage is \(O(N)\).
Footnotes and References
- Project Euler problem page: https://projecteuler.net/problem=677
- Combinatorial species: Wikipedia - Combinatorial species
- Cycle index: Wikipedia - Cycle index
- Burnside's lemma: Wikipedia - Burnside's lemma
- Dissymmetry theorem: Wikipedia - Dissymmetry theorem
Problem 677 source code
C++
#include <algorithm>
#include <cstdint>
#include <iostream>
#include <thread>
#include <vector>
using namespace std;
namespace {
constexpr int MOD = 1'000'000'007;
int 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 static_cast<int>(r);
}
int mod_norm(long long x) {
x %= MOD;
if (x < 0) x += MOD;
return static_cast<int>(x);
}
void update_square(const vector<int>& s, vector<int>& s_sq, int n, int limit_n) {
const int s_n = s[n];
if (s_n == 0) return;
const int max_k = min(n - 1, limit_n - n);
for (int k = 1; k <= max_k; ++k) {
const int v = s[k];
if (v == 0) continue;
const int idx = n + k;
const int add = static_cast<int>((2LL * s_n * v) % MOD);
int val = s_sq[idx] + add;
if (val >= MOD) val -= MOD;
s_sq[idx] = val;
}
if (2 * n <= limit_n) {
const int add = static_cast<int>((1LL * s_n * s_n) % MOD);
int val = s_sq[2 * n] + add;
if (val >= MOD) val -= MOD;
s_sq[2 * n] = val;
}
}
int conv_cube_at(const vector<int>& s, const vector<int>& s_sq, int t) {
if (t <= 1) return 0;
long long sum = 0;
for (int k = 1; k <= t - 1; ++k) {
const int a = s[k];
const int b = s_sq[t - k];
if (a == 0 || b == 0) continue;
sum += (1LL * a * b) % MOD;
}
return static_cast<int>(sum % MOD);
}
int conv_sc2_at(const vector<int>& s, int t) {
if (t <= 1) return 0;
long long sum = 0;
for (int u = 2; u <= t - 1; u += 2) {
const int k = t - u;
if (k < 1) break;
const int a = s[k];
const int b = s[u / 2];
if (a == 0 || b == 0) continue;
sum += (1LL * a * b) % MOD;
}
return static_cast<int>(sum % MOD);
}
vector<int> build_shift(const vector<int>& s, int k, int nmax) {
vector<int> out(nmax + 1, 0);
for (int i = 1; i * k <= nmax; ++i) {
out[i * k] = s[i];
}
return out;
}
vector<int> convolution_single(const vector<int>& a, const vector<int>& b, int nmax) {
vector<int> out(nmax + 1, 0);
for (int i = 1; i <= nmax; ++i) {
if (a[i] == 0) continue;
const int ai = a[i];
for (int j = 1; j + i <= nmax; ++j) {
if (b[j] == 0) continue;
const int idx = i + j;
const int add = static_cast<int>((1LL * ai * b[j]) % MOD);
int val = out[idx] + add;
if (val >= MOD) val -= MOD;
out[idx] = val;
}
}
return out;
}
vector<int> convolution_parallel(const vector<int>& a, const vector<int>& b, int nmax, unsigned threads) {
if (threads <= 1 || nmax < 512) {
return convolution_single(a, b, nmax);
}
vector<vector<int>> partials(threads, vector<int>(nmax + 1, 0));
vector<thread> workers;
const int chunk = (nmax + static_cast<int>(threads) - 1) / static_cast<int>(threads);
for (unsigned tid = 0; tid < threads; ++tid) {
int start = static_cast<int>(tid) * chunk + 1;
int end = min(nmax, start + chunk - 1);
workers.emplace_back([&, tid, start, end]() {
auto& local = partials[tid];
for (int i = start; i <= end; ++i) {
if (a[i] == 0) continue;
const int ai = a[i];
for (int j = 1; j + i <= nmax; ++j) {
if (b[j] == 0) continue;
const int idx = i + j;
const int add = static_cast<int>((1LL * ai * b[j]) % MOD);
int val = local[idx] + add;
if (val >= MOD) val -= MOD;
local[idx] = val;
}
}
});
}
for (auto& th : workers) th.join();
vector<int> out(nmax + 1, 0);
for (unsigned tid = 0; tid < threads; ++tid) {
const auto& local = partials[tid];
for (int i = 1; i <= nmax; ++i) {
int val = out[i] + local[i];
if (val >= MOD) val -= MOD;
out[i] = val;
}
}
return out;
}
vector<int> convolution_stride(const vector<int>& a, const vector<int>& base, int stride, int nmax, unsigned threads) {
const int max_k = nmax / stride;
if (threads <= 1 || max_k < 256) {
vector<int> out(nmax + 1, 0);
for (int k = 1; k <= max_k; ++k) {
const int coeff = base[k];
if (coeff == 0) continue;
const int offset = stride * k;
for (int j = 1; j + offset <= nmax; ++j) {
const int aj = a[j];
if (aj == 0) continue;
const int idx = j + offset;
const int add = static_cast<int>((1LL * aj * coeff) % MOD);
int val = out[idx] + add;
if (val >= MOD) val -= MOD;
out[idx] = val;
}
}
return out;
}
vector<vector<int>> partials(threads, vector<int>(nmax + 1, 0));
vector<thread> workers;
const int chunk = (max_k + static_cast<int>(threads) - 1) / static_cast<int>(threads);
for (unsigned tid = 0; tid < threads; ++tid) {
int start = static_cast<int>(tid) * chunk + 1;
int end = min(max_k, start + chunk - 1);
workers.emplace_back([&, tid, start, end]() {
auto& local = partials[tid];
for (int k = start; k <= end; ++k) {
const int coeff = base[k];
if (coeff == 0) continue;
const int offset = stride * k;
for (int j = 1; j + offset <= nmax; ++j) {
const int aj = a[j];
if (aj == 0) continue;
const int idx = j + offset;
const int add = static_cast<int>((1LL * aj * coeff) % MOD);
int val = local[idx] + add;
if (val >= MOD) val -= MOD;
local[idx] = val;
}
}
});
}
for (auto& th : workers) th.join();
vector<int> out(nmax + 1, 0);
for (unsigned tid = 0; tid < threads; ++tid) {
const auto& local = partials[tid];
for (int i = 1; i <= nmax; ++i) {
int val = out[i] + local[i];
if (val >= MOD) val -= MOD;
out[i] = val;
}
}
return out;
}
struct PlantedData {
vector<int> p_r;
vector<int> p_b;
vector<int> p_y;
vector<int> s_all;
vector<int> s_no_y;
vector<int> s_all_sq;
vector<int> s_no_y_sq;
};
PlantedData build_planted(int nmax) {
const int inv2 = mod_pow(2, MOD - 2);
const int inv6 = mod_pow(6, MOD - 2);
vector<int> p_r(nmax + 1, 0);
vector<int> p_b(nmax + 1, 0);
vector<int> p_y(nmax + 1, 0);
vector<int> s_all(nmax + 1, 0);
vector<int> s_no_y(nmax + 1, 0);
vector<int> s_all_sq(nmax + 1, 0);
vector<int> s_no_y_sq(nmax + 1, 0);
for (int n = 1; n <= nmax; ++n) {
const int t = n - 1;
const int base = (t == 0) ? 1 : 0;
const int s1_all = s_all[t];
const int s2_all = s_all_sq[t];
const int s3_all = conv_cube_at(s_all, s_all_sq, t);
const int sc2_all = conv_sc2_at(s_all, t);
const int c2_all = (t % 2 == 0) ? s_all[t / 2] : 0;
const int c3_all = (t % 3 == 0) ? s_all[t / 3] : 0;
const int m1_all = s1_all;
const int m2_all = static_cast<int>((1LL * (s2_all + c2_all) % MOD) * inv2 % MOD);
const int m3_all = static_cast<int>((1LL * (s3_all + 3LL * sc2_all + 2LL * c3_all) % MOD) * inv6 % MOD);
const int s1_no = s_no_y[t];
const int s2_no = s_no_y_sq[t];
const int s3_no = conv_cube_at(s_no_y, s_no_y_sq, t);
const int sc2_no = conv_sc2_at(s_no_y, t);
const int c2_no = (t % 2 == 0) ? s_no_y[t / 2] : 0;
const int c3_no = (t % 3 == 0) ? s_no_y[t / 3] : 0;
const int m1_no = s1_no;
const int m2_no = static_cast<int>((1LL * (s2_no + c2_no) % MOD) * inv2 % MOD);
const int m3_no = static_cast<int>((1LL * (s3_no + 3LL * sc2_no + 2LL * c3_no) % MOD) * inv6 % MOD);
p_r[n] = mod_norm(static_cast<long long>(base) + m1_all + m2_all + m3_all);
p_b[n] = mod_norm(static_cast<long long>(base) + m1_all + m2_all);
p_y[n] = mod_norm(static_cast<long long>(base) + m1_no + m2_no);
s_all[n] = mod_norm(static_cast<long long>(p_r[n]) + p_b[n] + p_y[n]);
s_no_y[n] = mod_norm(static_cast<long long>(p_r[n]) + p_b[n]);
update_square(s_all, s_all_sq, n, nmax);
update_square(s_no_y, s_no_y_sq, n, nmax);
}
PlantedData out;
out.p_r = std::move(p_r);
out.p_b = std::move(p_b);
out.p_y = std::move(p_y);
out.s_all = std::move(s_all);
out.s_no_y = std::move(s_no_y);
out.s_all_sq = std::move(s_all_sq);
out.s_no_y_sq = std::move(s_no_y_sq);
return out;
}
} // namespace
int main() {
const int nmax = 10000;
const int inv2 = mod_pow(2, MOD - 2);
const int inv6 = mod_pow(6, MOD - 2);
const int inv24 = mod_pow(24, MOD - 2);
unsigned threads = thread::hardware_concurrency();
if (threads == 0) threads = 1;
PlantedData planted = build_planted(nmax);
const vector<int>& p_all = planted.s_all;
const vector<int>& p_no_y = planted.s_no_y;
const vector<int>& p_y = planted.p_y;
vector<int> c2_all = build_shift(p_all, 2, nmax);
vector<int> c3_all = build_shift(p_all, 3, nmax);
vector<int> c4_all = build_shift(p_all, 4, nmax);
vector<int> c2_no = build_shift(p_no_y, 2, nmax);
vector<int> c3_no = build_shift(p_no_y, 3, nmax);
vector<int> all_cube = convolution_parallel(planted.s_all_sq, p_all, nmax, threads);
vector<int> all_four = convolution_parallel(planted.s_all_sq, planted.s_all_sq, nmax, threads);
vector<int> all_c2 = convolution_stride(p_all, p_all, 2, nmax, threads);
vector<int> all_sq_c2 = convolution_stride(planted.s_all_sq, p_all, 2, nmax, threads);
vector<int> all_c3 = convolution_stride(p_all, p_all, 3, nmax, threads);
vector<int> no_cube = convolution_parallel(planted.s_no_y_sq, p_no_y, nmax, threads);
vector<int> no_c2 = convolution_stride(p_no_y, p_no_y, 2, nmax, threads);
vector<int> c2_sq_all(nmax + 1, 0);
for (int i = 1; 2 * i <= nmax; ++i) {
c2_sq_all[2 * i] = planted.s_all_sq[i];
}
vector<int> m1_all = p_all;
vector<int> m2_all(nmax + 1, 0);
vector<int> m3_all(nmax + 1, 0);
vector<int> m4_all(nmax + 1, 0);
vector<int> m1_no = p_no_y;
vector<int> m2_no(nmax + 1, 0);
vector<int> m3_no(nmax + 1, 0);
for (int t = 0; t <= nmax; ++t) {
m2_all[t] = static_cast<int>((1LL * (planted.s_all_sq[t] + c2_all[t]) % MOD) * inv2 % MOD);
m2_no[t] = static_cast<int>((1LL * (planted.s_no_y_sq[t] + c2_no[t]) % MOD) * inv2 % MOD);
m3_all[t] = static_cast<int>((1LL * (all_cube[t] + 3LL * all_c2[t] + 2LL * c3_all[t]) % MOD) * inv6 % MOD);
m3_no[t] = static_cast<int>((1LL * (no_cube[t] + 3LL * no_c2[t] + 2LL * c3_no[t]) % MOD) * inv6 % MOD);
m4_all[t] = static_cast<int>((1LL * (all_four[t] +
6LL * all_sq_c2[t] +
3LL * c2_sq_all[t] +
8LL * all_c3[t] +
6LL * c4_all[t]) % MOD) * inv24 % MOD);
}
vector<int> a_all(nmax + 1, 0);
vector<int> a_r(nmax + 1, 0);
vector<int> a_b(nmax + 1, 0);
vector<int> a_y(nmax + 1, 0);
for (int n = 1; n <= nmax; ++n) {
const int t = n - 1;
const int base = (t == 0) ? 1 : 0;
a_r[n] = mod_norm(static_cast<long long>(base) + m1_all[t] + m2_all[t] + m3_all[t] + m4_all[t]);
a_b[n] = mod_norm(static_cast<long long>(base) + m1_all[t] + m2_all[t] + m3_all[t]);
a_y[n] = mod_norm(static_cast<long long>(base) + m1_no[t] + m2_no[t] + m3_no[t]);
a_all[n] = mod_norm(static_cast<long long>(a_r[n]) + a_b[n] + a_y[n]);
}
vector<int> p_y_sq = convolution_parallel(p_y, p_y, nmax, threads);
vector<int> p_all_x2 = build_shift(p_all, 2, nmax);
vector<int> p_y_x2 = build_shift(p_y, 2, nmax);
vector<int> g(nmax + 1, 0);
for (int n = 1; n <= nmax; ++n) {
const int d = mod_norm(static_cast<long long>(planted.s_all_sq[n]) - p_y_sq[n]);
const int e = mod_norm(static_cast<long long>(planted.s_all_sq[n]) + p_all_x2[n] - p_y_sq[n] - p_y_x2[n]);
const int e_half = static_cast<int>((1LL * e) * inv2 % MOD);
g[n] = mod_norm(static_cast<long long>(a_all[n]) + e_half - d);
}
struct Check {
int n;
int expected;
};
const Check checks[] = {
{2, 5},
{3, 15},
{4, 57},
{10, 710249},
{100, 919747298},
};
for (const auto& chk : checks) {
if (g[chk.n] != chk.expected) {
cerr << "Validation failed for g(" << chk.n << "): got "
<< g[chk.n] << " expected " << chk.expected << "\n";
return 1;
}
}
cout << g[nmax] << "\n";
return 0;
}
Python
def solve():
MOD = 10**9+7; NMAX = 10000
inv2 = pow(2,MOD-2,MOD); inv6 = pow(6,MOD-2,MOD); inv24 = pow(24,MOD-2,MOD)
# Build planted tree sequences
s_all=[0]*(NMAX+1); s_no=[0]*(NMAX+1)
s_all_sq=[0]*(NMAX+1); s_no_sq=[0]*(NMAX+1)
p_r=[0]*(NMAX+1); p_b=[0]*(NMAX+1); p_y=[0]*(NMAX+1)
def up_sq(s,sq,n):
v=s[n]
if v==0: return
mx=min(n-1,NMAX-n)
for k in range(1,mx+1):
if s[k]==0: continue
sq[n+k]=(sq[n+k]+2*v*s[k])%MOD
if 2*n<=NMAX: sq[2*n]=(sq[2*n]+v*v)%MOD
def c3(s,sq,t):
if t<=1: return 0
r=0
for k in range(1,t): r=(r+s[k]*sq[t-k])%MOD
return r
def sc2(s,t):
if t<=1: return 0
r=0
for u in range(2,t,2):
k=t-u
if k<1: break
r=(r+s[k]*s[u//2])%MOD
return r
for n in range(1,NMAX+1):
t=n-1; base=1 if t==0 else 0
s1a=s_all[t]; s2a=s_all_sq[t]; s3a=c3(s_all,s_all_sq,t)
sc2a=sc2(s_all,t); c2a=s_all[t//2] if t%2==0 else 0
c3a=s_all[t//3] if t%3==0 else 0
m1a=s1a; m2a=(s2a+c2a)*inv2%MOD
m3a=(s3a+3*sc2a+2*c3a)%MOD*inv6%MOD
s1n=s_no[t]; s2n=s_no_sq[t]; s3n=c3(s_no,s_no_sq,t)
sc2n=sc2(s_no,t); c2n=s_no[t//2] if t%2==0 else 0
c3n=s_no[t//3] if t%3==0 else 0
m1n=s1n; m2n=(s2n+c2n)*inv2%MOD
p_r[n]=(base+m1a+m2a+m3a)%MOD
p_b[n]=(base+m1a+m2a)%MOD
p_y[n]=(base+m1n+m2n)%MOD
s_all[n]=(p_r[n]+p_b[n]+p_y[n])%MOD
s_no[n]=(p_r[n]+p_b[n])%MOD
up_sq(s_all,s_all_sq,n); up_sq(s_no,s_no_sq,n)
def shift(s,k):
out=[0]*(NMAX+1)
for i in range(1,NMAX//k+1): out[i*k]=s[i]
return out
def conv(a,b):
out=[0]*(NMAX+1)
for i in range(1,NMAX+1):
if a[i]==0: continue
for j in range(1,NMAX-i+1):
if b[j]==0: continue
out[i+j]=(out[i+j]+a[i]*b[j])%MOD
return out
def conv_stride(a,base,stride):
out=[0]*(NMAX+1)
for k in range(1,NMAX//stride+1):
c=base[k]
if c==0: continue
off=stride*k
for j in range(1,NMAX-off+1):
if a[j]==0: continue
out[j+off]=(out[j+off]+a[j]*c)%MOD
return out
c2_all=shift(s_all,2); c3_all=shift(s_all,3); c4_all=shift(s_all,4)
c2_no=shift(s_no,2); c3_no=shift(s_no,3)
all_cube=conv(s_all_sq,s_all); all_four=conv(s_all_sq,s_all_sq)
all_c2=conv_stride(s_all,s_all,2); all_sq_c2=conv_stride(s_all_sq,s_all,2)
all_c3=conv_stride(s_all,s_all,3)
no_cube=conv(s_no_sq,s_no); no_c2=conv_stride(s_no,s_no,2)
c2_sq_all=[0]*(NMAX+1)
for i in range(1,NMAX//2+1): c2_sq_all[2*i]=s_all_sq[i]
a_all=[0]*(NMAX+1)
for n in range(1,NMAX+1):
t=n-1; base=1 if t==0 else 0
m2a=(s_all_sq[t]+c2_all[t])*inv2%MOD
m3a=(all_cube[t]+3*all_c2[t]+2*c3_all[t])%MOD*inv6%MOD
m4a=(all_four[t]+6*all_sq_c2[t]+3*c2_sq_all[t]+8*all_c3[t]+6*c4_all[t])%MOD*inv24%MOD
m2n=(s_no_sq[t]+c2_no[t])*inv2%MOD
m3n=(no_cube[t]+3*no_c2[t]+2*c3_no[t])%MOD*inv6%MOD
ar=(base+s_all[t]+m2a+m3a+m4a)%MOD
ab=(base+s_all[t]+m2a+m3a)%MOD
ay=(base+s_no[t]+m2n+m3n)%MOD
a_all[n]=(ar+ab+ay)%MOD
py_sq=conv(p_y,p_y)
pa_x2=shift(s_all,2); py_x2=shift(p_y,2)
g=[0]*(NMAX+1)
for n in range(1,NMAX+1):
d=(s_all_sq[n]-py_sq[n])%MOD
e2=(s_all_sq[n]+pa_x2[n]-py_sq[n]-py_x2[n])%MOD
eh=e2*inv2%MOD
g[n]=(a_all[n]+eh-d)%MOD
return str(g[NMAX]%MOD)
if __name__=='__main__':
print(solve())
Java
public class Euler677 {
static final int MOD = 1000000007;
static int modPow(long a, long e) {
long r = 1;
a %= MOD;
while (e > 0) {
if ((e & 1L) != 0)
r = (r * a) % MOD;
a = (a * a) % MOD;
e >>= 1L;
}
return (int) r;
}
static int modNorm(long x) {
x %= MOD;
if (x < 0)
x += MOD;
return (int) x;
}
static void updateSquare(int[] s, int[] sSq, int n, int limitN) {
int sn = s[n];
if (sn == 0)
return;
int maxK = Math.min(n - 1, limitN - n);
long add = (2L * sn) % MOD;
for (int k = 1; k <= maxK; ++k) {
if (s[k] != 0) {
sSq[n + k] = (int) ((sSq[n + k] + add * s[k]) % MOD);
}
}
if (2 * n <= limitN) {
sSq[2 * n] = (int) ((sSq[2 * n] + (long) sn * sn) % MOD);
}
}
static int convCubeAt(int[] s, int[] sSq, int t) {
if (t <= 1)
return 0;
long sum = 0;
for (int k = 1; k < t; ++k) {
if (s[k] == 0 || sSq[t - k] == 0)
continue;
sum += (long) s[k] * sSq[t - k];
if (sum >= 8000000000000000000L)
sum %= MOD; // protect from long overflow
}
return (int) (sum % MOD);
}
static int convSc2At(int[] s, int t) {
if (t <= 1)
return 0;
long sum = 0;
for (int u = 2; u < t; u += 2) {
int k = t - u;
if (s[k] == 0 || s[u / 2] == 0)
continue;
sum += (long) s[k] * s[u / 2];
if (sum >= 8000000000000000000L)
sum %= MOD;
}
return (int) (sum % MOD);
}
static int[] buildShift(int[] s, int k, int nmax) {
int[] out = new int[nmax + 1];
for (int i = 1; i * k <= nmax; ++i) {
out[i * k] = s[i];
}
return out;
}
static int[] convolutionSingle(int[] a, int[] b, int nmax) {
int[] out = new int[nmax + 1];
for (int i = 1; i <= nmax; ++i) {
if (a[i] == 0)
continue;
long ai = a[i];
int limit = nmax - i;
for (int j = 1; j <= limit; ++j) {
if (b[j] == 0)
continue;
long val = out[i + j] + ai * b[j];
out[i + j] = (int) (val % MOD);
}
}
return out;
}
static int[] convolutionStride(int[] a, int[] base, int stride, int nmax) {
int[] out = new int[nmax + 1];
int maxK = nmax / stride;
for (int k = 1; k <= maxK; ++k) {
long coeff = base[k];
if (coeff == 0)
continue;
int offset = stride * k;
int limit = nmax - offset;
for (int j = 1; j <= limit; ++j) {
if (a[j] == 0)
continue;
long val = out[j + offset] + a[j] * coeff;
out[j + offset] = (int) (val % MOD);
}
}
return out;
}
static class PlantedData {
int[] p_r, p_b, p_y, s_all, s_no_y, s_all_sq, s_no_y_sq;
}
static PlantedData buildPlanted(int nmax) {
int inv2 = modPow(2, MOD - 2);
int inv6 = modPow(6, MOD - 2);
int[] p_r = new int[nmax + 1];
int[] p_b = new int[nmax + 1];
int[] p_y = new int[nmax + 1];
int[] s_all = new int[nmax + 1];
int[] s_no_y = new int[nmax + 1];
int[] s_all_sq = new int[nmax + 1];
int[] s_no_y_sq = new int[nmax + 1];
for (int n = 1; n <= nmax; ++n) {
int t = n - 1;
int base = (t == 0) ? 1 : 0;
int s1_all = s_all[t];
int s2_all = s_all_sq[t];
int s3_all = convCubeAt(s_all, s_all_sq, t);
int sc2_all = convSc2At(s_all, t);
int c2_all = (t % 2 == 0) ? s_all[t / 2] : 0;
int c3_all = (t % 3 == 0) ? s_all[t / 3] : 0;
int m1_all = s1_all;
int m2_all = (int) ((1L * (s2_all + c2_all) % MOD) * inv2 % MOD);
int m3_all = (int) ((1L * (s3_all + 3L * sc2_all + 2L * c3_all) % MOD) * inv6 % MOD);
int s1_no = s_no_y[t];
int s2_no = s_no_y_sq[t];
int s3_no = convCubeAt(s_no_y, s_no_y_sq, t);
int sc2_no = convSc2At(s_no_y, t);
int c2_no = (t % 2 == 0) ? s_no_y[t / 2] : 0;
int c3_no = (t % 3 == 0) ? s_no_y[t / 3] : 0;
int m1_no = s1_no;
int m2_no = (int) ((1L * (s2_no + c2_no) % MOD) * inv2 % MOD);
int m3_no = (int) ((1L * (s3_no + 3L * sc2_no + 2L * c3_no) % MOD) * inv6 % MOD);
p_r[n] = modNorm((long) base + m1_all + m2_all + m3_all);
p_b[n] = modNorm((long) base + m1_all + m2_all);
p_y[n] = modNorm((long) base + m1_no + m2_no);
s_all[n] = modNorm((long) p_r[n] + p_b[n] + p_y[n]);
s_no_y[n] = modNorm((long) p_r[n] + p_b[n]);
updateSquare(s_all, s_all_sq, n, nmax);
updateSquare(s_no_y, s_no_y_sq, n, nmax);
}
PlantedData pd = new PlantedData();
pd.p_r = p_r;
pd.p_b = p_b;
pd.p_y = p_y;
pd.s_all = s_all;
pd.s_no_y = s_no_y;
pd.s_all_sq = s_all_sq;
pd.s_no_y_sq = s_no_y_sq;
return pd;
}
public static String solve() {
int nmax = 10000;
int inv2 = modPow(2, MOD - 2);
int inv6 = modPow(6, MOD - 2);
int inv24 = modPow(24, MOD - 2);
PlantedData pd = buildPlanted(nmax);
int[] p_all = pd.s_all;
int[] p_no_y = pd.s_no_y;
int[] p_y = pd.p_y;
int[] c2_all = buildShift(p_all, 2, nmax);
int[] c3_all = buildShift(p_all, 3, nmax);
int[] c4_all = buildShift(p_all, 4, nmax);
int[] c2_no = buildShift(p_no_y, 2, nmax);
int[] c3_no = buildShift(p_no_y, 3, nmax);
int[] all_cube = convolutionSingle(pd.s_all_sq, p_all, nmax);
int[] all_four = convolutionSingle(pd.s_all_sq, pd.s_all_sq, nmax);
int[] all_c2 = convolutionStride(p_all, p_all, 2, nmax);
int[] all_sq_c2 = convolutionStride(pd.s_all_sq, p_all, 2, nmax);
int[] all_c3 = convolutionStride(p_all, p_all, 3, nmax);
int[] no_cube = convolutionSingle(pd.s_no_y_sq, p_no_y, nmax);
int[] no_c2 = convolutionStride(p_no_y, p_no_y, 2, nmax);
int[] c2_sq_all = new int[nmax + 1];
for (int i = 1; 2 * i <= nmax; ++i) {
c2_sq_all[2 * i] = pd.s_all_sq[i];
}
int[] m1_all = p_all;
int[] m2_all = new int[nmax + 1];
int[] m3_all = new int[nmax + 1];
int[] m4_all = new int[nmax + 1];
int[] m1_no = p_no_y;
int[] m2_no = new int[nmax + 1];
int[] m3_no = new int[nmax + 1];
for (int t = 0; t <= nmax; ++t) {
m2_all[t] = (int) ((1L * (pd.s_all_sq[t] + c2_all[t]) % MOD) * inv2 % MOD);
m2_no[t] = (int) ((1L * (pd.s_no_y_sq[t] + c2_no[t]) % MOD) * inv2 % MOD);
m3_all[t] = (int) ((1L * (all_cube[t] + 3L * all_c2[t] + 2L * c3_all[t]) % MOD) * inv6 % MOD);
m3_no[t] = (int) ((1L * (no_cube[t] + 3L * no_c2[t] + 2L * c3_no[t]) % MOD) * inv6 % MOD);
m4_all[t] = (int) ((1L
* (all_four[t] + 6L * all_sq_c2[t] + 3L * c2_sq_all[t] + 8L * all_c3[t] + 6L * c4_all[t]) % MOD)
* inv24 % MOD);
}
int[] a_all = new int[nmax + 1];
int[] a_r = new int[nmax + 1];
int[] a_b = new int[nmax + 1];
int[] a_y = new int[nmax + 1];
for (int n = 1; n <= nmax; ++n) {
int t = n - 1;
int base = (t == 0) ? 1 : 0;
a_r[n] = modNorm((long) base + m1_all[t] + m2_all[t] + m3_all[t] + m4_all[t]);
a_b[n] = modNorm((long) base + m1_all[t] + m2_all[t] + m3_all[t]);
a_y[n] = modNorm((long) base + m1_no[t] + m2_no[t] + m3_no[t]);
a_all[n] = modNorm((long) a_r[n] + a_b[n] + a_y[n]);
}
int[] p_y_sq = convolutionSingle(p_y, p_y, nmax);
int[] p_all_x2 = buildShift(p_all, 2, nmax);
int[] p_y_x2 = buildShift(p_y, 2, nmax);
int[] g = new int[nmax + 1];
for (int n = 1; n <= nmax; ++n) {
int d = modNorm((long) pd.s_all_sq[n] - p_y_sq[n]);
int e = modNorm((long) pd.s_all_sq[n] + p_all_x2[n] - p_y_sq[n] - p_y_x2[n]);
int e_half = (int) ((1L * e) * inv2 % MOD);
g[n] = modNorm((long) a_all[n] + e_half - d);
}
return Integer.toString(g[nmax]);
}
public static void main(String[] args) {
System.out.println(solve());
}
}