Detailed mathematical approach

Problem Summary

The target quantity is a summatory arithmetic function \(G(n)\), evaluated modulo \(10^9+7\), with the Project Euler input at \(n=10^{12}\). The native implementations show that the summand is multiplicative and that its local behavior depends on whether an odd prime is congruent to \(1\) or \(3\) modulo \(4\).

The relevant Dirichlet character is

$\chi_4(m)=\begin{cases} 0,& 2\mid m,\\ 1,& m\equiv 1\pmod 4,\\ -1,& m\equiv 3\pmod 4. \end{cases}$

Reading the prime and prime-power contributions from the implementations gives the multiplicative term

$a(1)=1,\qquad a(p^e)=p^{3e}\left(1-\frac{\chi_4(p)}{p^3}\right)=p^{3e}-\chi_4(p)p^{3e-3}.$

Therefore the problem is to compute

$G(n)=\sum_{m\le n} a(m)\pmod{10^9+7}.$

A direct loop to \(10^{12}\) is not viable, so the implementation uses a compressed prime-prefix sieve together with a DFS over admissible prime-power products.

Mathematical Approach

The entire method can be reconstructed from two facts visible in the implementations: a prime first appears with weight \(p^3-\chi_4(p)\), and every extra copy of the same prime multiplies the contribution by another factor \(p^3\).

Step 1: Identify the Local Prime Weight

For a prime \(p\), the first occurrence contributes

$w(p)=p^3-\chi_4(p).$

So the three residue classes behave as follows:

$w(2)=8,\qquad w(p)=p^3-1\ \text{if }p\equiv 1\pmod 4,\qquad w(p)=p^3+1\ \text{if }p\equiv 3\pmod 4.$

This is the arithmetic signature of the problem: primes \(1\bmod 4\) and \(3\bmod 4\) are treated differently, exactly through the character \(\chi_4\).

Step 2: Extend from Primes to Prime Powers

Once a branch already contains \(p\), increasing the exponent by \(1\) multiplies the current local contribution by \(p^3\). Therefore

$a(p^e)=w(p)\,p^{3(e-1)}=(p^3-\chi_4(p))p^{3(e-1)}=p^{3e}-\chi_4(p)p^{3e-3}.$

In particular, the factor for \(p=2\) is simply \(2^{3e}\), because \(\chi_4(2)=0\). For odd primes the correction term is exactly one unit of \(p^{3e-3}\), with sign determined by \(p\bmod 4\).

Step 3: Rebuild the Global Multiplicative Function

Independent prime branches multiply, so for

$n=\prod_{p^e\parallel n} p^e$

we get

$a(n)=\prod_{p^e\parallel n} a(p^e)=n^3\prod_{p\mid n}\left(1-\frac{\chi_4(p)}{p^3}\right).$

This is the closed form encoded by the DFS. An equivalent divisor-sum identity is

$a(n)=\sum_{d\mid n}\mu(d)\chi_4(d)\left(\frac{n}{d}\right)^3,$

because only squarefree divisors survive the Möbius factor. That identity is not used directly in the code, but it confirms the prime-power formula and makes the multiplicativity transparent.

Step 4: Separate the Prime Contribution

The summatory function can be decomposed into the contribution of \(1\), the contribution of primes, and the remaining composite branches. The prime part is

$\sum_{p\le x} a(p)=\sum_{p\le x}\bigl(p^3-\chi_4(p)\bigr).$

For that reason the implementation builds two prime-prefix tables:

$S_3(x)=\sum_{p\le x} p^3,\qquad C_4(x)=\sum_{p\le x}\chi_4(p).$

Then the prime contribution is their difference

$S(x)=S_3(x)-C_4(x)=\sum_{p\le x}\bigl(p^3-\chi_4(p)\bigr).$

Both tables start from easy full-integer prefixes and then sieve away composite contributions prime by prime. For cubes, the starting formula is

$\sum_{m=2}^{x} m^3=\left(\frac{x(x+1)}{2}\right)^2-1.$

For the mod-\(4\) character, the full-integer prefix is periodic:

$\sum_{m=2}^{x}\chi_4(m)=\begin{cases} 0,& x\equiv 1,2\pmod 4,\\ -1,& x\equiv 0,3\pmod 4. \end{cases}$

Step 5: Use Floor-Division Compression and DFS

The values \(\left\lfloor n/i\right\rfloor\) repeat heavily, so the implementation stores prime-prefix information only on the compressed domain

$\mathcal{D}(n)=\left\{1,2,\dots,\left\lfloor\frac{n}{\lfloor\sqrt n\rfloor}\right\rfloor\right\}\cup\left\{\left\lfloor\frac{n}{i}\right\rfloor:1\le i\le \lfloor\sqrt n\rfloor\right\}.$

This reduces the table size from \(O(n)\) to \(O(\sqrt n)\). After the prime-only prefixes are available on that domain, a DFS enumerates multiplicative branches in increasing prime order. Ordering the primes this way prevents double-counting. At each branch:

$\text{first copy of }p \Rightarrow p^3-\chi_4(p),\qquad \text{each extra copy of }p \Rightarrow \times p^3.$

The prime-prefix tables instantly add the cases where a branch is extended by one larger prime, while recursion handles products containing several larger primes.

Worked Example: \(G(10)\)

The local formula gives

$\begin{aligned} a(1)&=1,\\ a(2)&=2^3=8,\\ a(3)&=3^3\left(1+\frac{1}{3^3}\right)=28,\\ a(4)&=4^3=64,\\ a(5)&=5^3\left(1-\frac{1}{5^3}\right)=124,\\ a(6)&=a(2)a(3)=224,\\ a(7)&=7^3\left(1+\frac{1}{7^3}\right)=344,\\ a(8)&=8^3=512,\\ a(9)&=9^3\left(1+\frac{1}{3^3}\right)=756,\\ a(10)&=a(2)a(5)=992. \end{aligned}$

Summing these values yields

$G(10)=1+8+28+64+124+224+344+512+756+992=3053,$

which matches the checkpoint used by the native implementations.

How the Code Works

The C++ and Java implementations first generate all primes up to \(\lfloor\sqrt n\rfloor\). They then build two compressed prime-prefix tables: one for \(\sum_{p\le x} p^3\), and one for \(\sum_{p\le x}\chi_4(p)\). Each table starts from the corresponding prefix over all integers and applies a prime-sieving transform so that only prime contributions remain.

After subtracting those tables, the implementation has fast access to \(\sum_{p\le x}(p^3-\chi_4(p))\) for every needed compressed argument \(x\). A DFS then enumerates composite multiplicative branches in increasing prime order. When a new prime enters a branch, the factor is \(p^3-\chi_4(p)\); when the same prime is repeated, the factor gains another multiplier \(p^3\). The prefix tables let the implementation add all one-more-prime extensions in bulk, while recursion handles deeper products.

The final answer is

$G(n)=1+\sum_{p\le n}\bigl(p^3-\chi_4(p)\bigr)+\text{composite DFS contribution}\pmod{10^9+7}.$

The Python implementation is intentionally thin: it compiles and launches the native solver when needed, then returns the parsed numeric result.

Complexity Analysis

The compressed floor-division domain has size \(O(\sqrt n)\), so the memory usage is \(O(\sqrt n)\). The two prime-prefix transforms run on that compressed domain rather than on all integers up to \(n\); in standard Min_25-style analysis this is roughly \(O(n^{3/4}/\log n)\) work. The DFS over admissible prime-power products is much smaller than a direct scan to \(n\) and is practical for the required input size. Overall the method is decisively sublinear in \(n\) and is designed specifically for \(n=10^{12}\).

Footnotes and References

1. Project Euler problem page: https://projecteuler.net/problem=715
2. Dirichlet character: Wikipedia — Dirichlet character
3. Multiplicative function: Wikipedia — Multiplicative function
4. Möbius function: Wikipedia — Möbius function
5. Min_25 sieve overview: OI Wiki — Min_25 sieve