Detailed mathematical approach

Problem Summary

The problem defines an integer \(S(n)\) by summing over all \(t=1,2,\dots,n\). The term indexed by \(t\) contributes only when the divisibility condition \((n-t)! \mid t^{\,n-t}\) holds, and then its value is

$(-1)^{n-t}\frac{t^{\,n-t}}{(n-t)!}.$

The goal is not to compute one isolated \(S(n)\), but the enormous prefix sum

$A(N)=\sum_{n=1}^{N} S(n)$

for \(N=10^{18}\), modulo \(10^9+7\). A direct evaluation is hopeless, because even checking the divisibility condition term by term would already be far too slow. The successful approach turns that divisibility test into a simple description of which \(t\)-values are allowed, then reduces everything to a short sum of power sums.

Mathematical Approach

The implementations revolve around one re-indexing step, one valuation argument, and one fast method for evaluating \(\sum_{k=1}^{x} k^m\) at huge \(x\).

Fix the gap \(m=n-t\)

Write

$m=n-t,\qquad t=n-m.$

Then \(m\) ranges from \(0\) to \(n-1\), and the defining sum becomes

$S(n)=\sum_{m=0}^{n-1}\mathbf{1}_{\,m!\mid (n-m)^m}\,(-1)^m\frac{(n-m)^m}{m!}.$

For the prefix sum, it is better to regard \(m\) as the outer variable and \(t\) as the free value. Since \(n=t+m\), we obtain

$A(N)=\sum_{m\ge 0}\frac{(-1)^m}{m!}\sum_{\substack{t\ge 1\\ t+m\le N\\ m!\mid t^m}} t^m.$

Now the whole problem is: for each fixed \(m\), characterize the integers \(t\) for which \(m!\) divides \(t^m\).

The divisibility criterion is exactly a primorial criterion

Define

$P_m=\prod_{\substack{p\le m\\ p\text{ prime}}} p,$

with the empty product interpreted as \(1\), so \(P_0=P_1=1\). This is the primorial up to \(m\).

The crucial fact used by all three implementations is

$m!\mid t^m \quad\Longleftrightarrow\quad P_m\mid t.$

Why is this true? For every prime \(p\le m\), the factorial \(m!\) contains at least one factor \(p\). So if \(p\nmid t\), then \(v_p(t^m)=0\lt v_p(m!)\), and divisibility fails immediately. Therefore every prime \(p\le m\) must divide \(t\), which means \(P_m\mid t\).

Conversely, suppose \(P_m\mid t\). Then every prime \(p\le m\) divides \(t\), so \(v_p(t)\ge 1\) and hence \(v_p(t^m)\ge m\). On the other hand, Legendre's formula gives

$v_p(m!)=\sum_{j\ge 1}\left\lfloor \frac{m}{p^j}\right\rfloor \le m,$

so \(v_p(t^m)\ge v_p(m!)\) for every prime \(p\le m\). No larger prime matters, because it does not divide \(m!\). Thus \(m!\mid t^m\).

This is the decisive simplification: the complicated-looking factorial divisibility test collapses to the simple statement that \(t\) is a multiple of one explicit primorial.

Turn admissible \(t\) into a power sum

Once we know that the admissible values are exactly \(t=kP_m\), the inner sum becomes static:

$A(N)=\sum_{m\ge 0}\frac{(-1)^m}{m!}\sum_{kP_m+m\le N}(kP_m)^m.$

Pull \(P_m^m\) out of the inner sum:

$A(N)=\sum_{m\ge 0}(-1)^m\frac{P_m^m}{m!}\sum_{k=1}^{\left\lfloor\frac{N-m}{P_m}\right\rfloor} k^m.$

So for each \(m\) we only need the power sum

$F_m(x)=\sum_{k=1}^{x} k^m,\qquad x=\left\lfloor\frac{N-m}{P_m}\right\rfloor.$

The outer loop is short because primorials grow very fast. As soon as \(P_m>N-m\), the floor becomes \(0\), and no later \(m\) can contribute.

Worked example: evaluating \(S(10)\)

This small case shows how the primorial filter really selects the surviving terms. Here

$S(10)=\sum_{t=1}^{10}\mathbf{1}_{(10-t)!\mid t^{\,10-t}}(-1)^{10-t}\frac{t^{\,10-t}}{(10-t)!}.$

Write \(m=10-t\).

For \(m=4\), we have \(P_4=2\cdot 3=6\). Then \(t=6\) is a multiple of \(6\), so the term survives and contributes

$(-1)^4\frac{6^4}{4!}=\frac{1296}{24}=54.$

For \(m=3\), we have \(P_3=6\), but now \(t=7\) is not a multiple of \(6\), so the term vanishes. For \(m=2\), \(P_2=2\) and \(t=8\) is allowed, giving

$(-1)^2\frac{8^2}{2!}=32.$

The terms with \(m=1\) and \(m=0\) always survive, contributing \(-9\) and \(1\). All larger \(m\) fail the divisibility test. Therefore

$S(10)=54+32-9+1=78.$

This example captures the full logic of the large computation: each \(m\) contributes through multiples of \(P_m\), and nothing else matters.

Evaluate \(F_m(x)\) at astronomically large \(x\)

For fixed \(m\), the power sum \(F_m(x)\) is a polynomial in \(x\) of degree \(m+1\) by Faulhaber's theorem. Therefore it is completely determined by the \(m+2\) values

$F_m(0),F_m(1),\dots,F_m(m+1).$

If we set \(d=m+2\) and \(y_i=F_m(i)\), then for any \(x\)

$F_m(x)=\sum_{i=0}^{d-1} y_i \prod_{\substack{0\le j\le d-1\\ j\ne i}}\frac{x-j}{i-j}.$

Because the interpolation nodes are consecutive integers, the denominator simplifies to

$\prod_{\substack{0\le j\le d-1\\ j\ne i}}(i-j)=(-1)^{d-1-i} i!(d-1-i)!.$

The implementations therefore precompute factorials and inverse factorials modulo \(10^9+7\), then use prefix and suffix products of \((x-j)\) to evaluate the numerator of every Lagrange basis term in linear time.

This is especially efficient here because the contributing \(m\)-values are tiny compared with the modulus. For the target \(N=10^{18}\), only \(m=0,1,\dots,52\) contribute, so the interpolation tables never exceed \(d=54\) entries.

How the Code Works

Exact checkpoint phase

Each implementation contains a small exact routine that evaluates \(S(n)\) directly with arbitrary-precision integers. That exact routine is used only for validation: it confirms \(S(1)=1\), \(S(3)=-1\), the prefix \(\sum_{n=1}^{10}S(n)=43\), and then compares several larger prefixes against the fast method. This guards the derivation before the large-\(N\) computation is trusted.

Fast modular phase

The optimized solver walks upward through \(m\), maintaining three quantities: the actual primorial \(P_m\) for the stopping test, its residue modulo \(10^9+7\), and \(m!\) modulo \(10^9+7\). For each \(m\), it computes

$(-1)^m P_m^m (m!)^{-1}\pmod{10^9+7}$

and multiplies it by the modular value of \(F_m\!\left(\left\lfloor\frac{N-m}{P_m}\right\rfloor\right)\).

The power sum is not expanded symbolically. Instead, the code first tabulates the prefix values \(F_m(0),F_m(1),\dots,F_m(m+1)\), then performs one Lagrange evaluation modulo the prime modulus. The C++, Python, and Java implementations all follow this same mathematical pipeline; they only differ in language-specific integer types.

Complexity Analysis

The outer summation is extremely short. For \(N=10^{18}\), the primorial up to \(53\) already exceeds \(N-53\), so only \(53\) values of \(m\) contribute: \(m=0\) through \(m=52\).

For each fixed \(m\), building the sample table \(F_m(0),\dots,F_m(m+1)\), the factorial tables, and the prefix/suffix products all costs \(O(m)\) modular operations. Summing over all contributing \(m\), the total work is therefore \(O(m_{\max}^2)\) with \(m_{\max}=52\), and the memory usage is \(O(m_{\max})\). In concrete terms, the problem becomes small because the hard part is not the size of \(N\), but the very slow growth of the set of relevant exponents \(m\).

Footnotes and References

1. Problem page: https://projecteuler.net/problem=969
2. Primorial: Wikipedia - Primorial
3. p-adic valuation and Legendre's formula: Wikipedia - Legendre's formula
4. Power sums and Faulhaber's theorem: Wikipedia - Faulhaber's formula
5. Lagrange interpolation: Wikipedia - Lagrange polynomial