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

1. Problem page: https://projecteuler.net/problem=690
2. Lobster graph: Wikipedia - Lobster graph
3. Euler transform: Wikipedia - Euler transform
4. Partition function: Wikipedia - Partition (number theory)
5. Generating function: Wikipedia - Generating function
6. Dirichlet convolution: Wikipedia - Dirichlet convolution