Problem 690: Tom and Jerry
View on Project EulerProject Euler Problem 690 Solution
EulerSolve provides an optimized solution for Project Euler Problem 690, Tom and Jerry, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary The computation is naturally split into two stages. First we build a sequence \(L_n\) of lobster-tree counts through formal power-series algebra. Then we apply the Euler transform to \(L_n\) and obtain the target sequence \(T_n\). The requested answer is \(T_{2019}\bmod 10^9+7\), so the real task is to derive the generating functions in a form that can be truncated efficiently up to degree \(2019\). Mathematical Approach Let \(M=10^9+7\). Write $$L(x)=\sum_{n\ge 1} L_nx^n,\qquad T(x)=\sum_{n\ge 0} T_nx^n.$$ The C++, Python, and Java implementations all follow the same algebraic pipeline. Step 1: Start from the partition generating function The base series is the ordinary generating function for integer partitions: $$P(x)=\prod_{m\ge 1}\frac{1}{1-x^m}=\sum_{n\ge 0} p(n)x^n.$$ Its first coefficients are $$P(x)=1+x+2x^2+3x^3+5x^4+7x^5+11x^6+15x^7+\cdots.$$ The implementations also use the shifted series $$X(x)=xP(x)=x+x^2+2x^3+3x^4+5x^5+7x^6+\cdots.$$ This partition baseline is the raw supply of branch-multiset data from which the lobster-tree series is assembled....
Detailed mathematical approach
Problem Summary
The computation is naturally split into two stages. First we build a sequence \(L_n\) of lobster-tree counts through formal power-series algebra. Then we apply the Euler transform to \(L_n\) and obtain the target sequence \(T_n\). The requested answer is \(T_{2019}\bmod 10^9+7\), so the real task is to derive the generating functions in a form that can be truncated efficiently up to degree \(2019\).
Mathematical Approach
Let \(M=10^9+7\). Write
$$L(x)=\sum_{n\ge 1} L_nx^n,\qquad T(x)=\sum_{n\ge 0} T_nx^n.$$
The C++, Python, and Java implementations all follow the same algebraic pipeline.
Step 1: Start from the partition generating function
The base series is the ordinary generating function for integer partitions:
$$P(x)=\prod_{m\ge 1}\frac{1}{1-x^m}=\sum_{n\ge 0} p(n)x^n.$$
Its first coefficients are
$$P(x)=1+x+2x^2+3x^3+5x^4+7x^5+11x^6+15x^7+\cdots.$$
The implementations also use the shifted series
$$X(x)=xP(x)=x+x^2+2x^3+3x^4+5x^5+7x^6+\cdots.$$
This partition baseline is the raw supply of branch-multiset data from which the lobster-tree series is assembled.
Step 2: Form the two main rational contributions
Define two auxiliary series
$$A(x)=P(x)-\frac{1}{1-x},\qquad B(x)=P(x^2)-\frac{1}{1-x^2}.$$
From them the implementation constructs
$$U(x)=\frac{A(x)^2}{1-X(x)},$$
$$V(x)=\frac{B(x)\left(1+X(x)\right)}{1-x^2P(x^2)}.$$
These are the two large rational pieces that feed the lobster-tree count. The code averages the two cases, so the combined contribution later appears as \(\frac{1}{2}(U(x)+V(x))\).
Step 3: Add the shift and subtract the correction term
The averaged rational contribution is shifted by two vertices, so it enters \(L(x)\) as
$$\frac{x^2}{2}\left(U(x)+V(x)\right).$$
Then the simpler family \(X(x)\) is added back in. Finally a correction series is subtracted. In the implementations the correction coefficients \(r_n\) satisfy
$$r_0=1,\qquad r_n=r_{n-1}+r_{n-2}-r_{n-3}\quad (n\ge 1),$$
with the missing negative indices treated as \(0\). Therefore the correction generating function is
$$R(x)=\sum_{n\ge 0} r_nx^n=\frac{1}{1-x-x^2+x^3}=\frac{1}{(1-x)(1-x^2)}.$$
So the lobster-tree generating function used by the implementations is
$$L(x)=X(x)+\frac{x^2}{2}\left(U(x)+V(x)\right)-x^3R(x).$$
Substituting the definitions of \(U\), \(V\), and \(R\) gives the closed form
$$L(x)=xP(x)+\frac{x^2}{2}\left(\frac{\left(P(x)-\frac{1}{1-x}\right)^2}{1-xP(x)}+\frac{\left(P(x^2)-\frac{1}{1-x^2}\right)\left(1+xP(x)\right)}{1-x^2P(x^2)}\right)-\frac{x^3}{(1-x)(1-x^2)}.$$
Step 4: Expand the first lobster coefficients
Expanding the rational pieces shows
$$U(x)=x^4+5x^5+18x^6+53x^7+\cdots,$$
$$V(x)=x^4+x^5+4x^6+5x^7+\cdots,$$
$$\frac{x^3}{(1-x)(1-x^2)}=x^3+x^4+2x^5+2x^6+3x^7+3x^8+\cdots.$$
Combining these with \(X(x)\) gives
$$L(x)=x+x^2+x^3+2x^4+3x^5+6x^6+11x^7+23x^8+47x^9+105x^{10}+\cdots.$$
Hence
$$L_1,\dots,L_{10}=1,1,1,2,3,6,11,23,47,105,$$
which matches the checked prefix in the implementations.
Step 5: Apply the Euler transform
The target series is the Euler transform of the lobster counts:
$$T(x)=\prod_{d\ge 1}(1-x^d)^{-L_d}.$$
Define the divisor sum
$$c_m=\sum_{d\mid m} d\,L_d.$$
Taking a logarithmic derivative gives
$$\log T(x)=\sum_{d\ge 1} L_d\sum_{j\ge 1}\frac{x^{dj}}{j},$$
so
$$\frac{xT'(x)}{T(x)}=\sum_{m\ge 1} c_mx^m.$$
Comparing coefficients with \(T(x)=\sum_{n\ge 0}T_nx^n\) yields
$$T_0=1,\qquad nT_n=\sum_{k=1}^{n} c_kT_{n-k}.$$
Since the modulus is prime, division by \(n\) is performed with modular inverses:
$$T_n=\frac{1}{n}\sum_{k=1}^{n} c_kT_{n-k}\pmod{M}.$$
Worked Example
Using the first lobster counts \(L_1=1\), \(L_2=1\), \(L_3=1\), \(L_4=2\), we obtain
$$c_1=1,\qquad c_2=1+2=3,\qquad c_3=1+3=4,\qquad c_4=1+2+8=11.$$
Now compute the first transformed values:
$$T_1=\frac{c_1T_0}{1}=1,$$
$$T_2=\frac{c_1T_1+c_2T_0}{2}=\frac{1+3}{2}=2,$$
$$T_3=\frac{c_1T_2+c_2T_1+c_3T_0}{3}=\frac{2+3+4}{3}=3,$$
$$T_4=\frac{c_1T_3+c_2T_2+c_3T_1+c_4T_0}{4}=\frac{3+6+4+11}{4}=6.$$
Continuing in the same way gives \(T_7=37\), \(T_{10}=328\), and \(T_{20}=1416269\), exactly the values used as checkpoints by the implementation.
How the Code Works
The C++, Python, and Java implementations precompute modular inverses \(1^{-1},2^{-1},\dots,n^{-1}\) because the formulas divide by \(2\) and by \(n\). They then build the partition series \(P(x)\) with the standard dynamic program for partition numbers. Next they create the two rational series above, invert the denominators coefficient-by-coefficient as formal power series, and multiply series with truncated convolutions so only terms up to degree \(n\) are kept.
After assembling \(L_1,\dots,L_n\), the implementation computes every divisor sum \(c_m\) by adding the contribution of each divisor to all of its multiples. Finally it evaluates the Euler-transform recurrence from \(T_0=1\) up to the requested index. The C++ version optionally parallelizes the largest convolution loops, but all three languages implement the same mathematics and return the same coefficients.
Complexity Analysis
Let \(n\) be the requested index. Building the partition series costs \(O(n^2)\) time and \(O(n)\) memory. Each truncated convolution and each formal power-series inversion is also \(O(n^2)\), and only a constant number of such arrays are stored at once. The divisor-sum stage costs \(O(n\log n)\), while the final Euler-transform recurrence is \(O(n^2)\). Therefore the overall running time is \(O(n^2)\) and the memory usage is \(O(n)\).
Footnotes and References
- Problem page: https://projecteuler.net/problem=690
- Lobster graph: Wikipedia - Lobster graph
- Euler transform: Wikipedia - Euler transform
- Partition function: Wikipedia - Partition (number theory)
- Generating function: Wikipedia - Generating function
- Dirichlet convolution: Wikipedia - Dirichlet convolution
Problem 690 source code
C++
#include <algorithm>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <vector>
namespace {
using i64 = std::int64_t;
using u32 = std::uint32_t;
constexpr i64 kMod = 1'000'000'007LL;
constexpr int kDefaultN = 2019;
struct Options {
int n = kDefaultN;
bool run_checkpoints = true;
bool allow_multithreading = true;
unsigned requested_threads = 0U;
};
inline i64 add_mod(const i64 a, const i64 b) {
i64 out = a + b;
if (out >= kMod) {
out -= kMod;
}
return out;
}
inline i64 sub_mod(const i64 a, const i64 b) {
i64 out = a - b;
if (out < 0) {
out += kMod;
}
return out;
}
inline i64 mul_mod(const i64 a, const i64 b) {
return static_cast<i64>((static_cast<__int128>(a) * static_cast<__int128>(b)) % kMod);
}
i64 pow_mod(i64 base, i64 exp) {
i64 out = 1;
while (exp > 0) {
if ((exp & 1LL) != 0LL) {
out = mul_mod(out, base);
}
base = mul_mod(base, base);
exp >>= 1LL;
}
return out;
}
bool parse_u32_after_prefix(const std::string& arg, const char* prefix, u32& value) {
const std::string p(prefix);
if (arg.rfind(p, 0) != 0U) {
return false;
}
const std::string tail = arg.substr(p.size());
if (tail.empty()) {
return false;
}
std::uint64_t parsed = 0ULL;
for (const char c : tail) {
if (c < '0' || c > '9') {
return false;
}
parsed = parsed * 10ULL + static_cast<std::uint64_t>(c - '0');
if (parsed > static_cast<std::uint64_t>(std::numeric_limits<u32>::max())) {
return false;
}
}
value = static_cast<u32>(parsed);
return true;
}
bool parse_unsigned_after_prefix(const std::string& arg,
const char* prefix,
unsigned& value) {
u32 parsed = 0U;
if (!parse_u32_after_prefix(arg, prefix, parsed)) {
return false;
}
value = static_cast<unsigned>(parsed);
return true;
}
bool parse_arguments(const 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 (arg == "--single-thread") {
options.allow_multithreading = false;
continue;
}
u32 parsed_u32 = 0U;
if (parse_u32_after_prefix(arg, "--n=", parsed_u32)) {
if (parsed_u32 > static_cast<u32>(std::numeric_limits<int>::max())) {
std::cerr << "--n is too large.\n";
return false;
}
options.n = static_cast<int>(parsed_u32);
continue;
}
unsigned parsed_unsigned = 0U;
if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
options.requested_threads = parsed_unsigned;
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
if (options.n < 1) {
std::cerr << "--n must be at least 1.\n";
return false;
}
return true;
}
unsigned choose_thread_count(const bool allow_multithreading,
const unsigned requested_threads,
const int workload) {
if (!allow_multithreading || workload < 512) {
return 1U;
}
unsigned threads = requested_threads;
if (threads == 0U) {
threads = std::thread::hardware_concurrency();
if (threads == 0U) {
threads = 1U;
}
}
return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload)));
}
std::vector<i64> precompute_inverses(const int n) {
std::vector<i64> inv(static_cast<std::size_t>(n) + 1ULL, 0LL);
if (n >= 1) {
inv[1] = 1LL;
}
for (int i = 2; i <= n; ++i) {
const i64 quotient = kMod / static_cast<i64>(i);
const i64 remainder = kMod % static_cast<i64>(i);
inv[static_cast<std::size_t>(i)] =
sub_mod(0LL, mul_mod(quotient, inv[static_cast<std::size_t>(remainder)]));
}
return inv;
}
std::vector<i64> convolve_truncated(const std::vector<i64>& a,
const std::vector<i64>& b,
const int n,
const bool allow_multithreading,
const unsigned requested_threads) {
std::vector<i64> out(static_cast<std::size_t>(n) + 1ULL, 0LL);
if (a.empty() || b.empty() || n < 0) {
return out;
}
const int max_i = std::min<int>(n, static_cast<int>(a.size()) - 1);
const int max_j_total = std::min<int>(n, static_cast<int>(b.size()) - 1);
if (max_i < 0 || max_j_total < 0) {
return out;
}
const unsigned threads = choose_thread_count(allow_multithreading, requested_threads, max_i + 1);
if (threads == 1U) {
for (int i = 0; i <= max_i; ++i) {
const i64 ai = a[static_cast<std::size_t>(i)];
if (ai == 0LL) {
continue;
}
const int max_j = std::min(max_j_total, n - i);
for (int j = 0; j <= max_j; ++j) {
const i64 bj = b[static_cast<std::size_t>(j)];
if (bj == 0LL) {
continue;
}
const std::size_t idx = static_cast<std::size_t>(i + j);
out[idx] = add_mod(out[idx], mul_mod(ai, bj));
}
}
return out;
}
std::vector<std::vector<i64>> partial(
threads, std::vector<i64>(static_cast<std::size_t>(n) + 1ULL, 0LL));
std::vector<std::thread> pool;
pool.reserve(threads);
for (unsigned t = 0; t < threads; ++t) {
const int start = static_cast<int>((static_cast<std::int64_t>(t) * (max_i + 1)) / threads);
const int end =
static_cast<int>((static_cast<std::int64_t>(t + 1U) * (max_i + 1)) / threads);
pool.emplace_back([&, t, start, end]() {
std::vector<i64>& local = partial[static_cast<std::size_t>(t)];
for (int i = start; i < end; ++i) {
const i64 ai = a[static_cast<std::size_t>(i)];
if (ai == 0LL) {
continue;
}
const int max_j = std::min(max_j_total, n - i);
for (int j = 0; j <= max_j; ++j) {
const i64 bj = b[static_cast<std::size_t>(j)];
if (bj == 0LL) {
continue;
}
const std::size_t idx = static_cast<std::size_t>(i + j);
local[idx] = add_mod(local[idx], mul_mod(ai, bj));
}
}
});
}
for (auto& th : pool) {
th.join();
}
for (unsigned t = 0; t < threads; ++t) {
const std::vector<i64>& local = partial[static_cast<std::size_t>(t)];
for (int i = 0; i <= n; ++i) {
out[static_cast<std::size_t>(i)] =
add_mod(out[static_cast<std::size_t>(i)], local[static_cast<std::size_t>(i)]);
}
}
return out;
}
std::vector<i64> invert_series(const std::vector<i64>& den, const int n) {
std::vector<i64> out(static_cast<std::size_t>(n) + 1ULL, 0LL);
out[0] = pow_mod(den[0], kMod - 2);
for (int i = 1; i <= n; ++i) {
i64 sum = 0LL;
for (int k = 1; k <= i; ++k) {
sum = add_mod(sum, mul_mod(den[static_cast<std::size_t>(k)],
out[static_cast<std::size_t>(i - k)]));
}
out[static_cast<std::size_t>(i)] = sub_mod(0LL, mul_mod(out[0], sum));
}
return out;
}
std::vector<i64> partition_series(const int n) {
std::vector<i64> p(static_cast<std::size_t>(n) + 1ULL, 0LL);
p[0] = 1LL;
for (int part = 1; part <= n; ++part) {
for (int i = part; i <= n; ++i) {
p[static_cast<std::size_t>(i)] =
add_mod(p[static_cast<std::size_t>(i)], p[static_cast<std::size_t>(i - part)]);
}
}
return p;
}
std::vector<i64> compute_lobster_tree_counts(const int n,
const std::vector<i64>& inv,
const bool allow_multithreading,
const unsigned requested_threads) {
const std::vector<i64> p = partition_series(n);
std::vector<i64> inv_1mx(static_cast<std::size_t>(n) + 1ULL, 1LL);
std::vector<i64> a(static_cast<std::size_t>(n) + 1ULL, 0LL);
for (int i = 0; i <= n; ++i) {
a[static_cast<std::size_t>(i)] =
sub_mod(p[static_cast<std::size_t>(i)], inv_1mx[static_cast<std::size_t>(i)]);
}
std::vector<i64> x_p(static_cast<std::size_t>(n) + 1ULL, 0LL);
for (int i = 1; i <= n; ++i) {
x_p[static_cast<std::size_t>(i)] = p[static_cast<std::size_t>(i - 1)];
}
std::vector<i64> den1(static_cast<std::size_t>(n) + 1ULL, 0LL);
den1[0] = 1LL;
for (int i = 1; i <= n; ++i) {
den1[static_cast<std::size_t>(i)] = sub_mod(0LL, x_p[static_cast<std::size_t>(i)]);
}
const std::vector<i64> inv_den1 = invert_series(den1, n);
const std::vector<i64> a_sq =
convolve_truncated(a, a, n, allow_multithreading, requested_threads);
const std::vector<i64> part1 =
convolve_truncated(a_sq, inv_den1, n, allow_multithreading, requested_threads);
std::vector<i64> p2(static_cast<std::size_t>(n) + 1ULL, 0LL);
for (int i = 0; i <= n; i += 2) {
p2[static_cast<std::size_t>(i)] = p[static_cast<std::size_t>(i / 2)];
}
std::vector<i64> inv_1mx2(static_cast<std::size_t>(n) + 1ULL, 0LL);
for (int i = 0; i <= n; i += 2) {
inv_1mx2[static_cast<std::size_t>(i)] = 1LL;
}
std::vector<i64> b(static_cast<std::size_t>(n) + 1ULL, 0LL);
for (int i = 0; i <= n; ++i) {
b[static_cast<std::size_t>(i)] =
sub_mod(p2[static_cast<std::size_t>(i)], inv_1mx2[static_cast<std::size_t>(i)]);
}
std::vector<i64> one_plus_x_p = x_p;
one_plus_x_p[0] = add_mod(one_plus_x_p[0], 1LL);
std::vector<i64> den2(static_cast<std::size_t>(n) + 1ULL, 0LL);
den2[0] = 1LL;
for (int i = 2; i <= n; ++i) {
den2[static_cast<std::size_t>(i)] = sub_mod(0LL, p2[static_cast<std::size_t>(i - 2)]);
}
const std::vector<i64> inv_den2 = invert_series(den2, n);
const std::vector<i64> tmp =
convolve_truncated(b, one_plus_x_p, n, allow_multithreading, requested_threads);
const std::vector<i64> part2 =
convolve_truncated(tmp, inv_den2, n, allow_multithreading, requested_threads);
std::vector<i64> lobsters(static_cast<std::size_t>(n) + 1ULL, 0LL);
for (int i = 0; i + 2 <= n; ++i) {
const i64 sum = add_mod(part1[static_cast<std::size_t>(i)], part2[static_cast<std::size_t>(i)]);
lobsters[static_cast<std::size_t>(i + 2)] =
add_mod(lobsters[static_cast<std::size_t>(i + 2)],
mul_mod(sum, inv[2]));
}
for (int i = 1; i <= n; ++i) {
lobsters[static_cast<std::size_t>(i)] =
add_mod(lobsters[static_cast<std::size_t>(i)], x_p[static_cast<std::size_t>(i)]);
}
std::vector<i64> rec(static_cast<std::size_t>(n) + 1ULL, 0LL);
rec[0] = 1LL;
for (int i = 1; i <= n; ++i) {
i64 v = rec[static_cast<std::size_t>(i - 1)];
if (i >= 2) {
v = add_mod(v, rec[static_cast<std::size_t>(i - 2)]);
}
if (i >= 3) {
v = sub_mod(v, rec[static_cast<std::size_t>(i - 3)]);
}
rec[static_cast<std::size_t>(i)] = v;
}
for (int i = 0; i + 3 <= n; ++i) {
lobsters[static_cast<std::size_t>(i + 3)] =
sub_mod(lobsters[static_cast<std::size_t>(i + 3)], rec[static_cast<std::size_t>(i)]);
}
return lobsters;
}
std::vector<i64> euler_transform_from_tree_counts(const std::vector<i64>& tree_counts,
const int n,
const std::vector<i64>& inv) {
std::vector<i64> c(static_cast<std::size_t>(n) + 1ULL, 0LL);
for (int d = 1; d <= n; ++d) {
const i64 val = mul_mod(static_cast<i64>(d), tree_counts[static_cast<std::size_t>(d)]);
for (int m = d; m <= n; m += d) {
c[static_cast<std::size_t>(m)] = add_mod(c[static_cast<std::size_t>(m)], val);
}
}
std::vector<i64> out(static_cast<std::size_t>(n) + 1ULL, 0LL);
out[0] = 1LL;
for (int i = 1; i <= n; ++i) {
i64 sum = 0LL;
for (int k = 1; k <= i; ++k) {
sum = add_mod(sum,
mul_mod(c[static_cast<std::size_t>(k)],
out[static_cast<std::size_t>(i - k)]));
}
out[static_cast<std::size_t>(i)] = mul_mod(sum, inv[static_cast<std::size_t>(i)]);
}
return out;
}
bool run_checkpoints() {
constexpr int kCheckN = 20;
const std::vector<i64> inv = precompute_inverses(kCheckN);
const std::vector<i64> lobster_trees =
compute_lobster_tree_counts(kCheckN, inv, false, 1U);
const std::vector<i64> expected_tree_prefix = {
0LL, 1LL, 1LL, 1LL, 2LL, 3LL, 6LL, 11LL, 23LL, 47LL, 105LL, 231LL, 532LL};
for (std::size_t i = 1; i < expected_tree_prefix.size(); ++i) {
if (lobster_trees[i] != expected_tree_prefix[i]) {
std::cerr << "Checkpoint failed: lobster tree count mismatch at n=" << i << ".\n";
return false;
}
}
const std::vector<i64> tom = euler_transform_from_tree_counts(lobster_trees, kCheckN, inv);
if (tom[3] != 3LL) {
std::cerr << "Checkpoint failed: T(3) != 3.\n";
return false;
}
if (tom[7] != 37LL) {
std::cerr << "Checkpoint failed: T(7) != 37.\n";
return false;
}
if (tom[10] != 328LL) {
std::cerr << "Checkpoint failed: T(10) != 328.\n";
return false;
}
if (tom[20] != 1'416'269LL) {
std::cerr << "Checkpoint failed: T(20) != 1416269.\n";
return false;
}
return true;
}
i64 solve(const int n, const bool allow_multithreading, const unsigned requested_threads) {
const std::vector<i64> inv = precompute_inverses(n);
const std::vector<i64> lobster_trees =
compute_lobster_tree_counts(n, inv, allow_multithreading, requested_threads);
const std::vector<i64> tom = euler_transform_from_tree_counts(lobster_trees, n, inv);
return tom[static_cast<std::size_t>(n)];
}
} // 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;
}
const i64 answer = solve(options.n, options.allow_multithreading, options.requested_threads);
std::cout << answer << '\n';
return 0;
}
Python
MOD = 1000000007
def add_mod(a, b):
return (a + b) % MOD
def sub_mod(a, b):
return (a - b + MOD) % MOD
def mul_mod(a, b):
return (a * b) % MOD
def pow_mod(base, exp):
return pow(base, exp, MOD)
def precompute_inverses(n):
inv = [0] * (n + 1)
if n >= 1:
inv[1] = 1
for i in range(2, n + 1):
inv[i] = MOD - mul_mod(MOD // i, inv[MOD % i])
return inv
def convolve_truncated(a, b, n):
out = [0] * (n + 1)
for i in range(min(n + 1, len(a))):
ai = a[i]
if ai == 0:
continue
for j in range(min(n - i + 1, len(b))):
bj = b[j]
if bj == 0:
continue
out[i + j] = (out[i + j] + ai * bj) % MOD
return out
def invert_series(den, n):
out = [0] * (n + 1)
out[0] = pow_mod(den[0], MOD - 2)
for i in range(1, n + 1):
sum_val = 0
for k in range(1, i + 1):
sum_val = (sum_val + den[k] * out[i - k]) % MOD
out[i] = (MOD - mul_mod(out[0], sum_val)) % MOD
return out
def partition_series(n):
p = [0] * (n + 1)
p[0] = 1
for part in range(1, n + 1):
for i in range(part, n + 1):
p[i] = (p[i] + p[i - part]) % MOD
return p
def compute_lobster_tree_counts(n, inv):
p = partition_series(n)
inv_1mx = [1] * (n + 1)
a = [(p[i] - inv_1mx[i]) % MOD for i in range(n + 1)]
x_p = [0] * (n + 1)
for i in range(1, n + 1):
x_p[i] = p[i - 1]
den1 = [0] * (n + 1)
den1[0] = 1
for i in range(1, n + 1):
den1[i] = (-x_p[i]) % MOD
inv_den1 = invert_series(den1, n)
a_sq = convolve_truncated(a, a, n)
part1 = convolve_truncated(a_sq, inv_den1, n)
p2 = [0] * (n + 1)
for i in range(0, n + 1, 2):
p2[i] = p[i // 2]
inv_1mx2 = [0] * (n + 1)
for i in range(0, n + 1, 2):
inv_1mx2[i] = 1
b = [(p2[i] - inv_1mx2[i]) % MOD for i in range(n + 1)]
one_plus_x_p = list(x_p)
one_plus_x_p[0] = (one_plus_x_p[0] + 1) % MOD
den2 = [0] * (n + 1)
den2[0] = 1
for i in range(2, n + 1):
den2[i] = (-p2[i - 2]) % MOD
inv_den2 = invert_series(den2, n)
tmp = convolve_truncated(b, one_plus_x_p, n)
part2 = convolve_truncated(tmp, inv_den2, n)
lobsters = [0] * (n + 1)
for i in range(n - 1):
sum_val = (part1[i] + part2[i]) % MOD
lobsters[i + 2] = (lobsters[i + 2] + sum_val * inv[2]) % MOD
for i in range(1, n + 1):
lobsters[i] = (lobsters[i] + x_p[i]) % MOD
rec = [0] * (n + 1)
rec[0] = 1
for i in range(1, n + 1):
v = rec[i - 1]
if i >= 2: v = (v + rec[i - 2]) % MOD
if i >= 3: v = (v - rec[i - 3]) % MOD
rec[i] = v % MOD
for i in range(n - 2):
lobsters[i + 3] = (lobsters[i + 3] - rec[i]) % MOD
return [(v % MOD + MOD) % MOD for v in lobsters]
def euler_transform(tree_counts, n, inv):
c = [0] * (n + 1)
for d in range(1, n + 1):
val = (d * tree_counts[d]) % MOD
for m in range(d, n + 1, d):
c[m] = (c[m] + val) % MOD
out = [0] * (n + 1)
out[0] = 1
for i in range(1, n + 1):
sum_val = 0
for k in range(1, i + 1):
sum_val = (sum_val + c[k] * out[i - k]) % MOD
out[i] = (sum_val * inv[i]) % MOD
return out
def solve():
n = 2019
inv = precompute_inverses(n)
lobsters = compute_lobster_tree_counts(n, inv)
tom = euler_transform(lobsters, n, inv)
return str(tom[n])
if __name__ == '__main__':
print(solve())
Java
public class Euler690 {
static final long MOD = 1000000007L;
static long powMod(long base, long exp) {
long out = 1;
base %= MOD;
while (exp > 0) {
if ((exp & 1) != 0)
out = (out * base) % MOD;
base = (base * base) % MOD;
exp >>= 1;
}
return out;
}
static long[] precomputeInverses(int n) {
long[] inv = new long[n + 1];
if (n >= 1)
inv[1] = 1;
for (int i = 2; i <= n; ++i) {
long quotient = MOD / i;
long remainder = MOD % i;
inv[i] = (MOD - (quotient * inv[(int) remainder]) % MOD) % MOD;
}
return inv;
}
static long[] convolveTruncated(long[] a, long[] b, int n) {
long[] out = new long[n + 1];
int maxI = Math.min(n, a.length - 1);
int maxJTotal = Math.min(n, b.length - 1);
for (int i = 0; i <= maxI; ++i) {
long ai = a[i];
if (ai == 0)
continue;
int maxJ = Math.min(maxJTotal, n - i);
for (int j = 0; j <= maxJ; ++j) {
long bj = b[j];
if (bj == 0)
continue;
out[i + j] = (out[i + j] + ai * bj) % MOD;
}
}
return out;
}
static long[] invertSeries(long[] den, int n) {
long[] out = new long[n + 1];
out[0] = powMod(den[0], MOD - 2);
for (int i = 1; i <= n; ++i) {
long sum = 0;
for (int k = 1; k <= i; ++k) {
sum = (sum + den[k] * out[i - k]) % MOD;
}
out[i] = (MOD - (out[0] * sum) % MOD) % MOD;
}
return out;
}
static long[] partitionSeries(int n) {
long[] p = new long[n + 1];
p[0] = 1;
for (int part = 1; part <= n; ++part) {
for (int i = part; i <= n; ++i) {
p[i] = (p[i] + p[i - part]) % MOD;
}
}
return p;
}
static long[] computeLobsterTreeCounts(int n, long[] inv) {
long[] p = partitionSeries(n);
long[] inv1mx = new long[n + 1];
for (int i = 0; i <= n; i++)
inv1mx[i] = 1;
long[] a = new long[n + 1];
for (int i = 0; i <= n; ++i) {
a[i] = (p[i] - inv1mx[i] + MOD) % MOD;
}
long[] xP = new long[n + 1];
for (int i = 1; i <= n; ++i) {
xP[i] = p[i - 1];
}
long[] den1 = new long[n + 1];
den1[0] = 1;
for (int i = 1; i <= n; ++i) {
den1[i] = (MOD - xP[i]) % MOD;
}
long[] invDen1 = invertSeries(den1, n);
long[] aSq = convolveTruncated(a, a, n);
long[] part1 = convolveTruncated(aSq, invDen1, n);
long[] p2 = new long[n + 1];
for (int i = 0; i <= n; i += 2) {
p2[i] = p[i / 2];
}
long[] inv1mx2 = new long[n + 1];
for (int i = 0; i <= n; i += 2) {
inv1mx2[i] = 1;
}
long[] b = new long[n + 1];
for (int i = 0; i <= n; ++i) {
b[i] = (p2[i] - inv1mx2[i] + MOD) % MOD;
}
long[] onePlusXP = xP.clone();
onePlusXP[0] = (onePlusXP[0] + 1) % MOD;
long[] den2 = new long[n + 1];
den2[0] = 1;
for (int i = 2; i <= n; ++i) {
den2[i] = (MOD - p2[i - 2]) % MOD;
}
long[] invDen2 = invertSeries(den2, n);
long[] tmp = convolveTruncated(b, onePlusXP, n);
long[] part2 = convolveTruncated(tmp, invDen2, n);
long[] lobsters = new long[n + 1];
for (int i = 0; i + 2 <= n; ++i) {
long sum = (part1[i] + part2[i]) % MOD;
lobsters[i + 2] = (lobsters[i + 2] + sum * inv[2]) % MOD;
}
for (int i = 1; i <= n; ++i) {
lobsters[i] = (lobsters[i] + xP[i]) % MOD;
}
long[] rec = new long[n + 1];
rec[0] = 1;
for (int i = 1; i <= n; ++i) {
long v = rec[i - 1];
if (i >= 2)
v = (v + rec[i - 2]) % MOD;
if (i >= 3)
v = (v - rec[i - 3] + MOD) % MOD;
rec[i] = v;
}
for (int i = 0; i + 3 <= n; ++i) {
lobsters[i + 3] = (lobsters[i + 3] - rec[i] + MOD) % MOD;
}
return lobsters;
}
static long[] eulerTransform(long[] treeCounts, int n, long[] inv) {
long[] c = new long[n + 1];
for (int d = 1; d <= n; ++d) {
long val = (d * treeCounts[d]) % MOD;
for (int m = d; m <= n; m += d) {
c[m] = (c[m] + val) % MOD;
}
}
long[] out = new long[n + 1];
out[0] = 1;
for (int i = 1; i <= n; ++i) {
long sum = 0;
for (int k = 1; k <= i; ++k) {
sum = (sum + c[k] * out[i - k]) % MOD;
}
out[i] = (sum * inv[i]) % MOD;
}
return out;
}
public static String solve() {
int n = 2019;
long[] inv = precomputeInverses(n);
long[] lobsters = computeLobsterTreeCounts(n, inv);
long[] tom = eulerTransform(lobsters, n, inv);
return Long.toString(tom[n]);
}
public static void main(String[] args) {
System.out.println(solve());
}
}