Problem 632: Square Prime Factors

Project Euler Problem 632 Solution

EulerSolve provides an optimized solution for Project Euler Problem 632, Square Prime Factors, with C++, Python, Java, and a step-by-step mathematical explanation.

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$....

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

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

Problem 632 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <cmath>
#include <iostream>
#include <numeric>
#include <vector>

using i64 = long long;
using u64 = unsigned long long;

static constexpr int MOD = 1'000'000'007;

static std::vector<int> primes_upto(int n) {
    std::vector<bool> is_comp(n / 2 + 1, false);
    std::vector<int> primes;
    primes.reserve(6'000'000);
    primes.push_back(2);
    for (int i = 3; (i64)i * i <= n; i += 2) {
        if (is_comp[i / 2]) continue;
        for (int j = i * i; j <= n; j += i * 2) is_comp[j / 2] = true;
    }
    for (int i = 3; i <= n; i += 2) {
        if (!is_comp[i / 2]) primes.push_back(i);
    }
    return primes;
}

static inline i64 isqrt_floor(i64 x) {
    i64 r = (i64)std::sqrt((long double)x);
    while ((r + 1) * (r + 1) <= x) ++r;
    while (r * r > x) --r;
    return r;
}

static std::array<u64, 10> compute_A(u64 N, const std::vector<uint8_t> &omega,
                                     const std::vector<uint8_t> &sqfree) {
    const int limit = (int)isqrt_floor((i64)N);
    std::array<u64, 10> A{};
    int i = 1;
    while (i <= limit) {
        const u64 q = N / (u64)i / (u64)i;
        int r = (int)isqrt_floor((i64)(N / q));
        if (r > limit) r = limit;
        while (r + 1 <= limit && N / (u64)(r + 1) / (u64)(r + 1) == q) ++r;
        while (r > limit || N / (u64)r / (u64)r != q) --r;

        for (int d = i; d <= r; ++d) {
            if (!sqfree[d]) continue;
            const int w = omega[d];
            if (w < 10) A[w] += q;
        }
        i = r + 1;
    }
    return A;
}

static std::array<u64, 10> compute_C(u64 N, const std::vector<uint8_t> &omega,
                                     const std::vector<uint8_t> &sqfree) {
    std::array<std::array<i64, 10>, 10> binom{};
    for (int n = 0; n < 10; ++n) {
        binom[n][0] = binom[n][n] = 1;
        for (int k = 1; k < n; ++k) binom[n][k] = binom[n - 1][k - 1] + binom[n - 1][k];
    }

    const auto A = compute_A(N, omega, sqfree);
    std::array<u64, 10> C{};
    for (int k = 0; k < 10; ++k) {
        __int128 s = 0;
        for (int r = k; r < 10; ++r) {
            const __int128 term = (__int128)A[r] * binom[r][k] * (((r - k) & 1) ? -1 : 1);
            s += term;
        }
        assert(s >= 0 && s <= (__int128)N);
        C[k] = (u64)s;
    }
    return C;
}

int main() {
    static constexpr int LIM = 100'000'000;
    const auto primes = primes_upto(LIM);

    std::vector<uint8_t> omega(LIM + 1, 0);
    std::vector<uint8_t> sqfree(LIM + 1, 1);
    sqfree[0] = 0;

    for (int p : primes) {
        for (int m = p; m <= LIM; m += p) ++omega[m];
        const i64 p2 = (i64)p * p;
        if (p2 <= LIM) {
            for (int m = (int)p2; m <= LIM; m += (int)p2) sqfree[m] = 0;
        }
    }

    {
        const auto C10 = compute_C(10ULL, omega, sqfree);
        assert(C10[0] == 7 && C10[1] == 3);
        for (int k = 2; k < 10; ++k) assert(C10[k] == 0);
    }
    {
        const auto C = compute_C(100ULL, omega, sqfree);
        assert(C[0] == 61 && C[1] == 36 && C[2] == 3);
    }
    {
        const auto C = compute_C(1'000ULL, omega, sqfree);
        assert(C[0] == 608 && C[1] == 343 && C[2] == 48 && C[3] == 1);
    }
    {
        const auto C = compute_C(10'000ULL, omega, sqfree);
        assert(C[0] == 6083 && C[1] == 3363 && C[2] == 533 && C[3] == 21);
    }
    {
        const auto C = compute_C(100'000ULL, omega, sqfree);
        assert(C[0] == 60794 && C[1] == 33562 && C[2] == 5345 && C[3] == 297 && C[4] == 2);
    }
    {
        const auto C = compute_C(1'000'000ULL, omega, sqfree);
        assert(C[0] == 607926 && C[1] == 335438 && C[2] == 53358 && C[3] == 3218 && C[4] == 60);
    }
    {
        const auto C = compute_C(10'000'000ULL, omega, sqfree);
        assert(C[0] == 6079291 && C[1] == 3353956 && C[2] == 533140 && C[3] == 32777 && C[4] == 834 &&
               C[5] == 2);
    }
    {
        const auto C = compute_C(100'000'000ULL, omega, sqfree);
        assert(C[0] == 60792694 && C[1] == 33539196 && C[2] == 5329747 && C[3] == 329028 && C[4] == 9257 &&
               C[5] == 78);
    }

    const u64 N = 10'000'000'000'000'000ULL;
    const auto C = compute_C(N, omega, sqfree);
    {
        __int128 s = 0;
        for (u64 x : C) s += x;
        assert(s == (__int128)N);
    }

    i64 ans = 1;
    for (u64 x : C) {
        if (x == 0) continue;
        ans = (i64)((__int128)ans * (x % MOD) % MOD);
    }
    std::cout << ans << "\n";
    return 0;
}

Python

import math

def solve():
    MOD = 1000000007
    LIM = 100000000
    N = 10000000000000000

    # Sieve primes
    half = LIM // 2 + 1
    comp = bytearray(half)
    primes = [2]
    for i in range(3, int(LIM**0.5)+1, 2):
        if not comp[i//2]:
            for j in range(i*i, LIM+1, 2*i): comp[j//2] = 1
    for i in range(3, LIM+1, 2):
        if not comp[i//2]: primes.append(i)

    omega = bytearray(LIM+1)
    sqfree = bytearray(b'\x01'*(LIM+1)); sqfree[0] = 0
    for p in primes:
        for m in range(p, LIM+1, p):
            if omega[m] < 255: omega[m] += 1
        p2 = p*p
        if p2 <= LIM:
            for m in range(p2, LIM+1, p2): sqfree[m] = 0

    def isqrt(x):
        r = int(math.isqrt(x))
        while (r+1)*(r+1) <= x: r += 1
        while r*r > x: r -= 1
        return r

    def compute_A():
        limit = isqrt(N)
        A = [0]*10; i = 1
        while i <= limit:
            q = N // (i*i)
            r = isqrt(N // q)
            if r > limit: r = limit
            while r+1 <= limit and N // ((r+1)*(r+1)) == q: r += 1
            while r > limit or N // (r*r) != q: r -= 1
            for d in range(i, r+1):
                if not sqfree[d]: continue
                w = omega[d]
                if w < 10: A[w] += q
            i = r + 1
        return A

    binom = [[0]*10 for _ in range(10)]
    for n in range(10):
        binom[n][0] = binom[n][n] = 1
        for k in range(1, n): binom[n][k] = binom[n-1][k-1] + binom[n-1][k]

    A = compute_A()
    C = [0]*10
    for k in range(10):
        s = 0
        for r in range(k, 10):
            s += A[r] * binom[r][k] * (1 if (r-k) % 2 == 0 else -1)
        C[k] = s

    ans = 1
    for x in C:
        if x == 0: continue
        ans = ans * (x % MOD) % MOD
    return str(ans)

if __name__ == '__main__':
    print(solve())

Java

import java.util.ArrayList;

public class Euler632 {
    static final int MOD = 1000000007;

    static ArrayList<Integer> primes_upto(int n) {
        ArrayList<Integer> primes = new ArrayList<>();
        if (n < 2)
            return primes;
        byte[] is_comp = new byte[n / 2 + 1];
        primes.add(2);
        for (int i = 3; (long) i * i <= n; i += 2) {
            if (is_comp[i / 2] == 0) {
                for (int j = i * i; j <= n; j += i * 2) {
                    is_comp[j / 2] = 1;
                }
            }
        }
        for (int i = 3; i <= n; i += 2) {
            if (is_comp[i / 2] == 0) {
                primes.add(i);
            }
        }
        return primes;
    }

    static long isqrt_floor(long x) {
        if (x < 0)
            return 0;
        long r = (long) Math.sqrt(x);
        while ((r + 1) * (r + 1) <= x)
            r++;
        while (r * r > x)
            r--;
        return r;
    }

    static long[] compute_A(long N, byte[] omega, byte[] sqfree) {
        int limit = (int) isqrt_floor(N);
        long[] A = new long[10];
        int i = 1;
        while (i <= limit) {
            long q = N / ((long) i * i);
            int r = (int) isqrt_floor(N / q);
            if (r > limit)
                r = limit;
            while (r + 1 <= limit && N / ((long) (r + 1) * (r + 1)) == q)
                r++;
            while (r > limit || N / ((long) r * r) != q)
                r--;

            for (int d = i; d <= r; ++d) {
                if (sqfree[d] == 0)
                    continue;
                int w = omega[d];
                if (w < 10)
                    A[w] += q;
            }
            i = r + 1;
        }
        return A;
    }

    static long[] compute_C(long N, byte[] omega, byte[] sqfree) {
        long[][] binom = new long[10][10];
        for (int n = 0; n < 10; ++n) {
            binom[n][0] = binom[n][n] = 1;
            for (int k = 1; k < n; ++k) {
                binom[n][k] = binom[n - 1][k - 1] + binom[n - 1][k];
            }
        }

        long[] A = compute_A(N, omega, sqfree);
        long[] C = new long[10];
        for (int k = 0; k < 10; ++k) {
            long s = 0;
            for (int r = k; r < 10; ++r) {
                long term = A[r] * binom[r][k] * (((r - k) % 2 != 0) ? -1 : 1);
                s += term;
            }
            C[k] = s;
        }
        return C;
    }

    public static String solve() {
        int LIM = 100000000;
        ArrayList<Integer> primes = primes_upto(LIM);
        byte[] omega = new byte[LIM + 1];
        byte[] sqfree = new byte[LIM + 1];
        for (int i = 1; i <= LIM; i++)
            sqfree[i] = 1;
        sqfree[0] = 0;

        for (int p : primes) {
            for (int m = p; m <= LIM; m += p)
                omega[m]++;
            long p2 = (long) p * p;
            if (p2 <= LIM) {
                for (int m = (int) p2; m <= LIM; m += (int) p2) {
                    sqfree[m] = 0;
                }
            }
        }

        long N = 10000000000000000L;
        long[] C = compute_C(N, omega, sqfree);

        long ans = 1;
        for (long x : C) {
            if (x == 0)
                continue;
            ans = (ans * (x % MOD)) % MOD;
        }
        return Long.toString(ans);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}