Detailed mathematical approach

Problem Summary

Problem 559 asks for the value of \(Q(50000)\) modulo the prime

$$M=1000000123.$$

The matrix-counting statement is encoded through auxiliary quantities \(P(k,r,n)\). The C++, Python, and Java implementations never try to enumerate the matrices directly. Instead, they transform the count into a coefficient-extraction problem built from inverse factorial powers, alternating signs, and one outer sum over \(k\).

Mathematical Approach

Fix \(n\), \(r\), and \(k\). The implementations first compute a normalized core value and multiply by \((n!)^r\) only at the end. That normalization turns the combinatorial problem into a clean reciprocal-series recurrence.

Step 1: Normalize by factorial powers

Because \(M\) is prime and the computation only needs factorials up to \(n=50000\lt M\), every \(m!\) is invertible modulo \(M\). Define

$$u_m=(m!)^{-r}\pmod{M}\qquad (0\le m\le n).$$

These numbers are the basic weights of the method. All divisions by factorials are replaced by modular inverses, and the missing factor \((n!)^r\) is restored only after the normalized coefficient has been computed.

Step 2: Split \(n\) into full \(k\)-blocks and a remainder

Write

$$n=bk+s,\qquad b=\left\lfloor\frac{n}{k}\right\rfloor,\qquad 0\le s\lt k.$$

The recurrence advances in jumps of size \(k\). The quotient \(b\) counts how many full \(k\)-steps fit into \(n\), while the remainder \(s\) records the final incomplete part. This is why the formulas separate naturally into a full-block part and a remainder correction.

Step 3: Build the reciprocal power series

Introduce the formal power series

$$F_k(z)=\sum_{j\ge 0}(-1)^j u_{jk}z^j=1-u_k z+u_{2k}z^2-u_{3k}z^3+\cdots.$$

Now define its reciprocal

$$A_k(z)=\frac{1}{F_k(z)}=\sum_{i\ge 0} a_i z^i.$$

Comparing coefficients in \(A_k(z)F_k(z)=1\) gives

$$a_0=1,$$

$$a_i=\sum_{j=1}^{i}(-1)^{j+1}a_{i-j}u_{jk}\qquad (i\ge 1).$$

This alternating convolution is the heart of the algorithm. Once the sequence \(a_0,a_1,\dots,a_b\) is known, the rest is just a short finishing sum.

Step 4: Correct the final partial block

If \(s=0\), the normalized core contribution is simply

$$\gamma(n,k,r)=a_b.$$

If \(s\gt 0\), one shifted series is needed:

$$G_{k,s}(z)=\sum_{j\ge 0}(-1)^j u_{s+jk}z^j.$$

The desired core is then the coefficient of \(z^b\) in the product \(A_k(z)G_{k,s}(z)\):

$$\gamma(n,k,r)=[z^b](A_k(z)G_{k,s}(z))=\sum_{j=0}^{b}(-1)^j a_{b-j}u_{s+jk}.$$

So the same recurrence handles all complete \(k\)-blocks, and the shifted tail series accounts for the leftover \(s\) entries.

Step 5: Recover \(P(k,r,n)\) and assemble \(Q(n)\)

After the normalized core has been computed, the required count is

$$P(k,r,n)=(n!)^r\gamma(n,k,r)\pmod{M}.$$

For the Project Euler target, the implementations set \(r=n\) and sum over every \(k\):

$$Q(n)=(n!)^n\sum_{k=1}^{n}\gamma(n,k,n)\pmod{M}.$$

The full problem is therefore reduced to evaluating one reciprocal-series coefficient for each \(k\), summing those normalized cores, and multiplying once at the end.

Worked Example: \(P(1,2,3)=19\)

Take \(k=1\), \(r=2\), and \(n=3\). Then \(b=3\) and \(s=0\). The relevant weights are

$$u_1=\frac{1}{(1!)^2}=1,\qquad u_2=\frac{1}{(2!)^2}=\frac14,\qquad u_3=\frac{1}{(3!)^2}=\frac1{36}.$$

Now apply the recurrence:

$$a_0=1,$$

$$a_1=a_0u_1=1,$$

$$a_2=a_1u_1-a_0u_2=1-\frac14=\frac34,$$

$$a_3=a_2u_1-a_1u_2+a_0u_3=\frac34-\frac14+\frac1{36}=\frac{19}{36}.$$

Because \(s=0\), we have \(\gamma(3,1,2)=a_3\). Therefore

$$P(1,2,3)=(3!)^2\cdot \frac{19}{36}=36\cdot \frac{19}{36}=19,$$

which matches the built-in checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations follow the same mathematics. They first precompute factorials and inverse factorials, then raise each inverse factorial to the exponent \(r\) to obtain the table \(u_m\). For a fixed \(k\), they fill the coefficient sequence \(a_0,a_1,\dots,a_b\) with the alternating convolution above and evaluate the shifted finishing sum when \(s>0\).

To compute \(P(k,r,n)\), the implementation multiplies the normalized core by \((n!)^r\). To compute \(Q(n)\), it uses \(r=n\), repeats the same per-\(k\) core computation for every \(1\le k\le n\), adds the cores, and only then multiplies by \((n!)^n\). The C++ version parallelizes the outer loop over \(k\), the Java version runs the same logic serially, and the Python version acts as a thin bridge that invokes the compiled C++ solver and returns its numeric result.

Complexity Analysis

Precomputing factorials and inverse factorials costs \(O(n)\). Building the table \(u_m=(m!)^{-r}\) with repeated fast exponentiation costs \(O(n\log r)\). For one fixed \(k\), the recurrence length is \(b=\lfloor n/k\rfloor\), and the nested alternating sum costs \(O(b^2)\) time. Therefore

$$\sum_{k=1}^{n} O\!\left(\left\lfloor\frac{n}{k}\right\rfloor^2\right)=O(n^2),$$

so \(Q(n)\) is evaluated in overall \(O(n^2)\) time. The main tables use \(O(n)\) memory; the multithreaded C++ version adds one reusable work array per worker thread.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=559
2. Permutation: Wikipedia — Permutation
3. Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
4. Formal power series: Wikipedia — Formal power series
5. Modular multiplicative inverse: Wikipedia — Modular multiplicative inverse