Detailed mathematical approach

Problem Summary

For each integer \(n\le N\), let \(t(n)\) denote the number of distinct primes whose squares divide \(n\). In other words, \(t(n)\) counts how many different primes \(p\) satisfy \(p^2\mid n\).

We then define

$C_k(N)=\left|\left\{n\in\mathbb{Z}_{>0}: n\le N,\ t(n)=k\right\}\right|.$

The task for Problem 632 is to evaluate, at \(N=10^{16}\), the product of all nonzero values \(C_k(N)\) modulo \(10^9+7\). Since scanning all integers up to \(10^{16}\) is impossible, the solution counts them indirectly by grouping together the squarefree products of the primes whose squares divide each integer.

Mathematical Approach

The key idea is to first count integers by how many chosen squared prime factors they contain, and only afterwards recover the exact buckets \(C_k(N)\).

Step 1: Encode a chosen set of squared primes by a squarefree number

If an integer \(n\) is divisible by the squares of the distinct primes \(p_1,\dots,p_r\), then the squarefree product

$d=p_1p_2\cdots p_r$

satisfies \(d^2\mid n\). Conversely, every squarefree \(d\) with \(d^2\mid n\) identifies a set of distinct squared prime divisors of \(n\). Let \(\omega(d)\) denote the number of distinct prime factors of \(d\), and let \(S_r(N)\) be the set of squarefree integers \(d\) such that \(d^2\le N\) and \(\omega(d)=r\).

Step 2: Define the aggregated counts \(A_r(N)\)

For each \(r\ge 0\), define

$A_r(N)=\sum_{d\in S_r(N)}\left\lfloor\frac{N}{d^2}\right\rfloor.$

Each term counts the integers \(n\le N\) divisible by \(d^2\). Therefore \(A_r(N)\) is not yet the number of integers with exactly \(r\) squared prime factors. Instead, it is an overcount: an integer is counted once for every squarefree choice of \(r\) squared primes contained in it.

Step 3: Relate the overcounts to the exact counts

Suppose \(t(n)=t\). Then \(n\) has exactly \(t\) distinct primes whose squares divide it. To contribute to \(A_r(N)\), we may choose any \(r\) of those \(t\) primes, so this single integer contributes exactly

$\binom{t}{r}$

times to \(A_r(N)\). Summing over all integers with the same value of \(t\) yields the triangular identity

$A_r(N)=\sum_{t\ge r}\binom{t}{r}C_t(N).$

This is the central combinatorial statement used by the implementations.

Step 4: Recover \(C_k(N)\) by binomial inversion

The previous system is lower triangular with diagonal entries equal to \(1\), so it can be inverted explicitly:

$C_k(N)=\sum_{r\ge k}(-1)^{r-k}\binom{r}{k}A_r(N).$

Once the values \(A_r(N)\) are known, each exact bucket \(C_k(N)\) follows from a short alternating sum.

Step 5: Why only finitely many buckets can be nonzero

If \(t(n)\ge 9\), then \(n\) must be divisible by the square of the product of the first nine primes:

$\left(2\cdot3\cdot5\cdot7\cdot11\cdot13\cdot17\cdot19\cdot23\right)^2>10^{16}.$

That is impossible when \(n\le 10^{16}\). Therefore only \(k=0,1,\dots,8\) can occur for the target input. This explains why the implementations only need a very small fixed number of counting buckets.

Step 6: Worked Example for \(N=100\)

This small case makes the counting mechanism transparent and matches the checkpoint values used by the implementation.

First, \(A_0(100)=100\), because the only squarefree \(d\) with \(\omega(d)=0\) is \(d=1\).

For \(r=1\), the possible values are \(d=2,3,5,7\), so

$A_1(100)=\left\lfloor\frac{100}{4}\right\rfloor+\left\lfloor\frac{100}{9}\right\rfloor+\left\lfloor\frac{100}{25}\right\rfloor+\left\lfloor\frac{100}{49}\right\rfloor=25+11+4+2=42.$

For \(r=2\), the only squarefree products with square at most \(100\) are \(d=6\) and \(d=10\), hence

$A_2(100)=\left\lfloor\frac{100}{36}\right\rfloor+\left\lfloor\frac{100}{100}\right\rfloor=2+1=3.$

All higher \(A_r(100)\) vanish. Now invert the system:

$C_2(100)=A_2(100)=3,$

$C_1(100)=A_1(100)-2A_2(100)=42-6=36,$

$C_0(100)=A_0(100)-A_1(100)+A_2(100)=100-42+3=61.$

So among the integers up to \(100\), exactly \(61\) have no squared prime factor, \(36\) have exactly one squared prime factor, and \(3\) have exactly two. The total \(61+36+3=100\) is a useful consistency check.

How the Code Works

The C++, Python, and Java implementations all begin by setting the limit to \(\lfloor\sqrt{N}\rfloor=10^8\), because every relevant squarefree \(d\) must satisfy \(d^2\le N\). They generate all primes up to that bound and then fill two compact arrays over \(1\le d\le 10^8\): one stores how many distinct prime divisors each \(d\) has, and the other records whether \(d\) is squarefree.

After that preprocessing phase, the implementation scans all \(d\le \sqrt{N}\). Whenever \(d\) is squarefree, it adds \(\left\lfloor N/d^2\right\rfloor\) to the bucket corresponding to \(\omega(d)\). Consecutive values of \(d\) that share the same quotient \(\left\lfloor N/d^2\right\rfloor\) are handled together, so the expensive division work is reused across whole blocks instead of repeated blindly for every bound adjustment.

Once the aggregated buckets \(A_r(N)\) are available, the implementation builds a tiny Pascal triangle of binomial coefficients and applies the inversion formula above to recover the exact values \(C_k(N)\). Finally, it multiplies all positive counts modulo \(10^9+7\). The C++ implementation also verifies several small checkpoints such as \(N=10,100,10^3,\dots,10^8\) and confirms that the recovered buckets sum back to \(N\).

Complexity Analysis

Let \(L=\lfloor\sqrt{N}\rfloor\). The sieve phase and the marking of multiples and square-multiples take \(O(L\log\log L)\) time in the standard sieve model and \(O(L)\) memory. The accumulation pass over \(d\le L\) is \(O(L)\), and the binomial inversion is constant-time because only a fixed number of buckets can be nonzero. For the target input \(N=10^{16}\), this means working up to \(L=10^8\): large, but still feasible for a one-shot optimized computation.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=632
2. Square-free integer: Wikipedia — Square-free integer
3. Binomial transform and inversion: Wikipedia — Binomial transform
4. Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
5. Prime omega function: Wikipedia — Prime omega function