Detailed mathematical approach

Problem Summary

For each positive integer \(k\), define

$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k,$

where \(\operatorname{rad}(n)\) is the product of the distinct prime divisors of \(n\). The problem asks for

$\sum_{k=1}^{50} S_k(10^{12}) \pmod{10^9+7}.$

The direct approach would require evaluating \(\operatorname{rad}(n)\) for every \(n\le 10^{12}\), which is hopeless. The implementation instead starts from the easier power sum \(\sum n^k\), then corrects it prime by prime until the local behavior matches \(\operatorname{rad}(n)^k\).

Mathematical Approach

Fix one value of \(k\). The task is to compute \(S_k(N)\), and then repeat that computation for \(k=1,2,\dots,50\).

Step 1: Start from the easier sum \(\sum n^k\)

The multiplicative function \(n^k\) has local values

$n^k=\prod_{p^a\parallel n} p^{ak}.$

By contrast,

$\operatorname{rad}(n)^k=\prod_{p^a\parallel n} p^k.$

So the only difference is what happens when a prime appears with exponent \(a\ge 2\). If a prime occurs only once, both functions contribute the same factor \(p^k\).

The implementation therefore begins with the summatory function

$P_k(x)=\sum_{n\le x} n^k,$

which can be evaluated quickly for many different \(x\) by a Faulhaber polynomial modulo \(10^9+7\).

Step 2: Correct one prime at a time

Suppose we have a temporary multiplicative weight in which the primes already processed behave like \(\operatorname{rad}(n)^k\), while the remaining primes still behave like \(n^k\). For the next prime \(p\), only prime powers \(p^a\) with \(a\ge 2\) need correction.

For one local factor we want to replace

$p^{ak}\quad\text{by}\quad p^k \qquad (a\ge 2).$

The difference is

$p^{ak}-p^k=p^k\left(p^{(a-1)k}-1\right).$

Using a geometric series, this becomes

$p^{ak}-p^k=p^k(p^k-1)\left(1+p^k+p^{2k}+\cdots+p^{(a-2)k}\right).$

Therefore the update coefficient for prime \(p\) is

$t_p=p^k(p^k-1).$

Step 3: Derive the summatory recurrence

Let \(F(x)\) be the current summatory function before correcting prime \(p\). Then subtracting

$t_p\sum_{a\ge 2} F\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right)$

exactly changes every local factor \(p^{ak}\) with \(a\ge 2\) into \(p^k\).

To see why, consider an integer \(n=p^e m\) with \(p\nmid m\). Its old local contribution at \(p\) is \(p^{ek}\). The subtraction contributes

$t_p\sum_{a=2}^{e} p^{(e-a)k}=p^k(p^k-1)\sum_{j=0}^{e-2} p^{jk}=p^{ek}-p^k.$

After subtraction, the new local contribution is therefore

$p^{ek}-(p^{ek}-p^k)=p^k,$

which is exactly what \(\operatorname{rad}(p^e)^k\) requires. If \(e=0\) or \(e=1\), no subtraction occurs, and the factor is already correct.

So the prime-by-prime transition is

$F_{\text{new}}(x)=F_{\text{old}}(x)-t_p\sum_{a\ge 2} F_{\text{old}}\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right).$

Step 4: Only primes up to \(\sqrt{N}\) matter

If \(p>\sqrt{N}\), then \(p^2>N\). Such a prime can appear in any \(n\le N\) only with exponent \(0\) or \(1\). But for exponents \(0\) and \(1\), the weights \(n^k\) and \(\operatorname{rad}(n)^k\) already agree. Hence corrections are needed only for primes

$p\le \sqrt{N}.$

This is why the implementation sieves primes only up to \(\lfloor\sqrt{N}\rfloor\).

Step 5: Compress all reachable arguments

The recurrence never asks for arbitrary values of \(x\). Starting from \(N\), every transition divides by some \(p^a\) with \(a\ge 2\). After several corrections, every queried argument has the form

$\left\lfloor\frac{N}{m}\right\rfloor,$

where \(m\) is a powerful number, meaning that every prime exponent in \(m\) is either \(0\) or at least \(2\).

This is the crucial compression step. Instead of storing values for all \(1\le x\le N\), the implementation stores values only for the distinct quotients produced by powerful denominators. These quotients are sorted once, indexed once, and then reused for all fifty values of \(k\).

Step 6: Sum over \(k=1,2,\dots,50\)

For each fixed \(k\), the algorithm starts from \(P_k(x)\), applies the prime corrections in increasing prime order, and obtains

$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k \pmod{10^9+7}.$

The final answer is simply

$\sum_{k=1}^{50} S_k(N)\pmod{10^9+7}.$

Worked Example: \(S_1(10)=41\)

For \(k=1\), the base sum is

$P_1(10)=1+2+\cdots+10=55.$

Only primes up to \(\sqrt{10}\) need correction, so only \(p=2\) and \(p=3\) matter.

For \(p=2\),

$t_2=2(2-1)=2.$

The correction uses \(2^2\) and \(2^3\):

$2\left(P_1\!\left(\left\lfloor\frac{10}{4}\right\rfloor\right)+P_1\!\left(\left\lfloor\frac{10}{8}\right\rfloor\right)\right)=2\left(P_1(2)+P_1(1)\right)=2(3+1)=8.$

So the total becomes \(55-8=47\).

For \(p=3\),

$t_3=3(3-1)=6,$

and only \(3^2\) contributes:

$6\,P_1\!\left(\left\lfloor\frac{10}{9}\right\rfloor\right)=6\,P_1(1)=6.$

Therefore

$S_1(10)=55-8-6=41.$

This matches the direct check

$\operatorname{rad}(1),\dots,\operatorname{rad}(10)=1,2,3,2,5,6,7,2,3,10,$

whose sum is \(41\).

How the Code Works

The C++, Python, and Java implementations follow the same structure. First they generate every prime up to \(\lfloor\sqrt{N}\rfloor\). Then they precompute Faulhaber-style polynomial coefficients modulo \(10^9+7\), so that \(\sum_{n\le x} n^k\) can be evaluated quickly for any needed quotient \(x\).

Next they enumerate the powerful denominators that can arise in the recurrence. From those denominators they build the distinct quotient states \(\left\lfloor N/m\right\rfloor\), sort them in descending order, and create index tables so every future division lands on an already known state.

For each \(k\in\{1,\dots,50\}\), every stored state is initialized with the corresponding power sum \(P_k(x)\). The implementation then processes primes in increasing order. For a given prime \(p\), it computes the factor \(t_p=p^k(p^k-1)\), and for every state \(x\ge p^2\) it subtracts the indexed states at \(\left\lfloor x/p^2\right\rfloor,\left\lfloor x/p^3\right\rfloor,\dots\).

The states are updated from larger quotients to smaller quotients. That ordering is important: when the implementation needs the value at \(\left\lfloor x/p^a\right\rfloor\), it still reads the pre-update value for the current prime, which is exactly what the recurrence requires.

After all relevant primes have been processed, the state corresponding to \(x=N\) is \(S_k(N)\). The implementation adds that value to the running total and continues with the next \(k\).

Complexity Analysis

Sieving primes up to \(\sqrt{N}\) costs \(O(\sqrt{N}\log\log N)\) time and \(O(\sqrt{N})\) memory. The Faulhaber preprocessing for \(k\le 50\) is tiny compared with the rest of the computation.

The main optimization is that the recurrence is evaluated only on compressed quotient states of the form \(\left\lfloor N/m\right\rfloor\) with powerful \(m\). That state set is far smaller than \(\{1,2,\dots,N\}\), which is what makes \(N=10^{12}\) feasible.

For each \(k\), the runtime is driven by the number of reachable pairs consisting of a quotient state and a prime power \(p^a\) with \(a\ge 2\). In other words, the algorithm scales with the compressed transition graph rather than with all integers up to \(N\). Memory is proportional to the number of stored quotient states together with the prime list and a few auxiliary tables.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=639
2. Radical of an integer: Wikipedia - Radical of an integer
3. Faulhaber's formula: Wikipedia - Faulhaber's formula
4. Bernoulli numbers: Wikipedia - Bernoulli number
5. Powerful numbers: Wikipedia - Powerful number
6. Multiplicative functions: Wikipedia - Multiplicative function