Problem 668: Square Root Smooth Numbers

View on Project Euler

Project Euler Problem 668 Solution

EulerSolve provides an optimized solution for Project Euler Problem 668, Square Root Smooth Numbers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For \(1\le x\le N\), the problem calls \(x\) square root smooth when every prime factor of \(x\) is strictly smaller than \(\sqrt{x}\). The number \(1\) is included as well, because it has no prime factors at all. The goal is to compute the count \(C(N)\) for \(N=10^{10}\) without testing every integer individually. Mathematical Approach Let \(P^+(x)\) denote the largest prime factor of \(x\) for \(x\ge 2\). Instead of counting square root smooth integers directly, it is easier to count the complement and subtract from \(N\). Step 1: Count the Complement Define $$NS(N)=\#\left\{x\in\{1,\dots,N\}:x\ge 2,\ P^+(x)\ge \sqrt{x}\right\}.$$ Then the required answer is simply $$C(N)=N-NS(N).$$ So the whole problem becomes: how many integers up to \(N\) have a largest prime factor at least as large as their square root? Step 2: Unique Representation \(x=pq\) Take any non-smooth integer \(x\), and let \(p=P^+(x)\) be its largest prime factor. Write $$x=pq.$$ Because \(p\ge \sqrt{x}\), we immediately get $$q=\frac{x}{p}\le p.$$ Also \(x\le N\) implies $$q\le \left\lfloor\frac{N}{p}\right\rfloor.$$ Conversely, suppose \(p\) is prime and $$1\le q\le \min\left(p,\left\lfloor\frac{N}{p}\right\rfloor\right).$$ Then every prime factor of \(q\) is at most \(q\le p\), so the largest prime factor of \(pq\) is exactly \(p\)....

Detailed mathematical approach

Problem Summary

For \(1\le x\le N\), the problem calls \(x\) square root smooth when every prime factor of \(x\) is strictly smaller than \(\sqrt{x}\). The number \(1\) is included as well, because it has no prime factors at all. The goal is to compute the count \(C(N)\) for \(N=10^{10}\) without testing every integer individually.

Mathematical Approach

Let \(P^+(x)\) denote the largest prime factor of \(x\) for \(x\ge 2\). Instead of counting square root smooth integers directly, it is easier to count the complement and subtract from \(N\).

Step 1: Count the Complement

Define

$$NS(N)=\#\left\{x\in\{1,\dots,N\}:x\ge 2,\ P^+(x)\ge \sqrt{x}\right\}.$$

Then the required answer is simply

$$C(N)=N-NS(N).$$

So the whole problem becomes: how many integers up to \(N\) have a largest prime factor at least as large as their square root?

Step 2: Unique Representation \(x=pq\)

Take any non-smooth integer \(x\), and let \(p=P^+(x)\) be its largest prime factor. Write

$$x=pq.$$

Because \(p\ge \sqrt{x}\), we immediately get

$$q=\frac{x}{p}\le p.$$

Also \(x\le N\) implies

$$q\le \left\lfloor\frac{N}{p}\right\rfloor.$$

Conversely, suppose \(p\) is prime and

$$1\le q\le \min\left(p,\left\lfloor\frac{N}{p}\right\rfloor\right).$$

Then every prime factor of \(q\) is at most \(q\le p\), so the largest prime factor of \(pq\) is exactly \(p\). Since \(q\le p\), we also have \(p\ge \sqrt{pq}\), so \(pq\) is non-smooth. Therefore every non-smooth number is counted exactly once by a pair \((p,q)\), and

$$NS(N)=\sum_{p\le N}\min\left(p,\left\lfloor\frac{N}{p}\right\rfloor\right),$$

where the sum runs over primes \(p\).

Step 3: Re-index by the Smaller Factor

The same pairs can be counted by fixing \(q\) first. Because \(q\le p\) and \(pq\le N\), necessarily

$$q\le \sqrt{N}.$$

For a fixed \(q\), the admissible primes satisfy

$$q\le p\le \left\lfloor\frac{N}{q}\right\rfloor.$$

Hence

$$NS(N)=\sum_{q=1}^{\lfloor\sqrt N\rfloor}\left(\pi\left(\left\lfloor\frac{N}{q}\right\rfloor\right)-\pi(q-1)\right).$$

This is the form used by the Python and Java implementations. Once prime counts are known on the relevant arguments, the final summation is immediate.

Step 4: Equivalent Split at \(\sqrt N\)

From the prime-based sum we can also separate the cases \(p\le \sqrt N\) and \(p>\sqrt N\):

$$NS(N)=\sum_{p\le \sqrt N}p+\sum_{\sqrt N<p\le N}\left\lfloor\frac{N}{p}\right\rfloor.$$

In the second sum the quotient

$$k=\left\lfloor\frac{N}{p}\right\rfloor$$

is constant on intervals of primes. The relevant block is

$$\max\left(\left\lfloor\sqrt N\right\rfloor+1,\left\lfloor\frac{N}{k+1}\right\rfloor+1\right)\le p\le \left\lfloor\frac{N}{k}\right\rfloor,$$

so its contribution is

$$k\left(\pi\left(\left\lfloor\frac{N}{k}\right\rfloor\right)-\pi\left(\max\left(\left\lfloor\sqrt N\right\rfloor,\left\lfloor\frac{N}{k+1}\right\rfloor\right)\right)\right).$$

This is the viewpoint used by the C++ implementation: it groups large primes by equal quotient instead of iterating over every prime one by one.

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

Here \(\lfloor\sqrt{100}\rfloor=10\). The small-prime part contributes

$$2+3+5+7=17.$$

The large-prime blocks are

$$\begin{aligned} k=1&:\quad 1\cdot(\pi(100)-\pi(50))=10,\\ k=2&:\quad 2\cdot(\pi(50)-\pi(33))=8,\\ k=3&:\quad 3\cdot(\pi(33)-\pi(25))=6,\\ k=4&:\quad 4\cdot(\pi(25)-\pi(20))=4,\\ k=5&:\quad 5\cdot(\pi(20)-\pi(16))=10,\\ k=7&:\quad 7\cdot(\pi(14)-\pi(12))=7,\\ k=9&:\quad 9\cdot(\pi(11)-\pi(10))=9. \end{aligned}$$

The missing values \(k=6\) and \(k=8\) contribute \(0\) because the corresponding intervals contain no primes. Therefore

$$NS(100)=17+10+8+6+4+10+7+9=71,$$

and thus

$$C(100)=100-71=29.$$

This matches the checkpoint used by the implementations.

Step 6: Why Prime Counting Solves the Problem

After the transformation above, the original smooth-number question no longer requires factoring every integer up to \(N\). Everything reduces to evaluating \(\pi(x)\) on roughly \(2\sqrt N\) relevant arguments and combining those values with one of the two equivalent summation formulas.

How the Code Works

The C++, Python, and Java implementations all begin with the complement identity \(C(N)=N-NS(N)\). They differ only in how they obtain the needed prime-counting values \(\pi(x)\).

The C++ implementation uses a standard sieve for small values and a cached Lehmer prime-counting method for large arguments. It evaluates the split formula: first add all primes up to \(\lfloor\sqrt N\rfloor\), then process the quotient blocks where \(\left\lfloor N/p\right\rfloor\) is constant and multiply each quotient by the number of primes in that interval.

The Python and Java implementations precompute \(\pi(v)\) simultaneously on the distinct values

$$v\in\left\{\left\lfloor\frac{N}{1}\right\rfloor,\left\lfloor\frac{N}{2}\right\rfloor,\dots,\left\lfloor\frac{N}{\lfloor\sqrt N\rfloor}\right\rfloor\right\}\cup\{1,2,\dots,\lfloor\sqrt N\rfloor\}.$$

They use the standard combinatorial prime-counting update over all primes up to \(\lfloor\sqrt N\rfloor\), turning an initial table \(v-1\) into exact values of \(\pi(v)\). With those values available, they evaluate

$$NS(N)=\sum_{q=1}^{\lfloor\sqrt N\rfloor}\left(\pi\left(\left\lfloor\frac{N}{q}\right\rfloor\right)-\pi(q-1)\right)$$

and subtract the result from \(N\).

Complexity Analysis

The mathematical reduction leaves only \(O(\sqrt N)\) relevant \(q\)-values or quotient blocks. The Python and Java implementations store prime-counting data on about \(2\sqrt N\) distinct arguments, so their memory usage is \(O(\sqrt N)\), and their preprocessing is the usual sublinear combinatorial prime-counting pass over those values. The C++ implementation also performs only \(O(\sqrt N)\) outer accumulation steps, with the main cost coming from memoized prime-counting queries. In all three cases the method is vastly faster than scanning every integer up to \(N\) and is easily practical for \(N=10^{10}\).

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=668
  2. Largest prime factor: Wikipedia — Largest prime factor
  3. Prime-counting function: Wikipedia — Prime-counting function
  4. Smooth number: Wikipedia — Smooth number
  5. Sieve of Eratosthenes: Wikipedia — Sieve of Eratosthenes

Problem 668 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <unordered_map>
#include <vector>

namespace {

using u64 = std::uint64_t;
using i64 = std::int64_t;
using u128 = __uint128_t;

u64 isqrt_u64(u64 n) {
    u64 x = static_cast<u64>(std::sqrt(static_cast<long double>(n)));
    while ((x + 1ULL) <= n / (x + 1ULL)) {
        ++x;
    }
    while (x > 0ULL && x > n / x) {
        --x;
    }
    return x;
}

u64 icbrt_u64(u64 n) {
    u64 x = static_cast<u64>(std::cbrt(static_cast<long double>(n)));
    while ((x + 1ULL) <= n / ((x + 1ULL) * (x + 1ULL))) {
        ++x;
    }
    while (x > 0ULL && x > n / (x * x)) {
        --x;
    }
    return x;
}

u64 iroot4_u64(u64 n) {
    u64 x = static_cast<u64>(
        std::sqrt(static_cast<long double>(std::sqrt(static_cast<long double>(n)))));
    auto pow4 = [](u64 y) -> u128 { return static_cast<u128>(y) * y * y * y; };
    while (pow4(x + 1ULL) <= n) {
        ++x;
    }
    while (x > 0ULL && pow4(x) > n) {
        --x;
    }
    return x;
}

class PrimeCounting {
public:
    PrimeCounting() { build_sieve(5'000'000); }

    u64 pi(u64 n) {
        if (n <= static_cast<u64>(sieve_limit_)) {
            return pi_small_[static_cast<std::size_t>(n)];
        }
        auto it = pi_cache_.find(n);
        if (it != pi_cache_.end()) {
            return it->second;
        }

        const u64 a = pi(iroot4_u64(n));
        const u64 b = pi(isqrt_u64(n));
        const u64 c = pi(icbrt_u64(n));

        i64 sum = static_cast<i64>(phi(n, static_cast<int>(a))) +
                  static_cast<i64>((b + a - 2ULL) * (b - a + 1ULL) / 2ULL);

        for (u64 i = a + 1ULL; i <= b; ++i) {
            const u64 p = static_cast<u64>(primes_[static_cast<std::size_t>(i - 1ULL)]);
            const u64 w = n / p;
            sum -= static_cast<i64>(pi(w));
            if (i <= c) {
                const u64 lim = pi(isqrt_u64(w));
                for (u64 j = i; j <= lim; ++j) {
                    const u64 pj = static_cast<u64>(primes_[static_cast<std::size_t>(j - 1ULL)]);
                    sum -= static_cast<i64>(pi(w / pj) - (j - 1ULL));
                }
            }
        }

        const u64 out = static_cast<u64>(sum);
        pi_cache_.emplace(n, out);
        return out;
    }

    const std::vector<int>& primes() const { return primes_; }

private:
    void build_sieve(int limit) {
        sieve_limit_ = limit;
        std::vector<bool> is_comp(static_cast<std::size_t>(limit + 1), false);
        pi_small_.assign(static_cast<std::size_t>(limit + 1), 0ULL);

        for (int i = 2; i <= limit; ++i) {
            if (!is_comp[static_cast<std::size_t>(i)]) {
                primes_.push_back(i);
                if (i <= limit / i) {
                    for (int j = i * i; j <= limit; j += i) {
                        is_comp[static_cast<std::size_t>(j)] = true;
                    }
                }
            }
            pi_small_[static_cast<std::size_t>(i)] =
                pi_small_[static_cast<std::size_t>(i - 1)] +
                (!is_comp[static_cast<std::size_t>(i)] ? 1ULL : 0ULL);
        }

        pi_cache_.reserve(1 << 20);
        phi_cache_.reserve(1 << 20);
    }

    u64 phi(u64 x, int s) {
        if (s == 0) {
            return x;
        }
        if (s == 1) {
            return x - x / 2ULL;
        }
        if (s == 2) {
            return x - x / 2ULL - x / 3ULL + x / 6ULL;
        }
        if (s == 3) {
            return x - x / 2ULL - x / 3ULL - x / 5ULL + x / 6ULL + x / 10ULL + x / 15ULL -
                   x / 30ULL;
        }
        if (x <= static_cast<u64>(sieve_limit_) &&
            static_cast<u64>(primes_[static_cast<std::size_t>(s - 1)]) >= x) {
            return 1ULL;
        }

        const u64 key = (x << 6U) ^ static_cast<u64>(s);
        auto it = phi_cache_.find(key);
        if (it != phi_cache_.end()) {
            return it->second;
        }

        const u64 p = static_cast<u64>(primes_[static_cast<std::size_t>(s - 1)]);
        const u64 out = phi(x, s - 1) - phi(x / p, s - 1);
        phi_cache_.emplace(key, out);
        return out;
    }

    int sieve_limit_ = 0;
    std::vector<int> primes_;
    std::vector<u64> pi_small_;
    std::unordered_map<u64, u64> pi_cache_;
    std::unordered_map<u64, u64> phi_cache_;
};

u64 count_smooth(u64 n, PrimeCounting& pc) {
    const u64 root = isqrt_u64(n);

    u128 non_smooth = 0;

    for (int p : pc.primes()) {
        if (static_cast<u64>(p) > root) {
            break;
        }
        non_smooth += static_cast<u64>(p);
    }

    const u64 q_max = n / (root + 1ULL);
    for (u64 q = 1; q <= q_max; ++q) {
        const u64 left = n / (q + 1ULL) + 1ULL;
        const u64 right = n / q;
        const u64 l = std::max<u64>(root + 1ULL, left);
        if (l > right) {
            continue;
        }
        const u64 cnt = pc.pi(right) - pc.pi(l - 1ULL);
        non_smooth += static_cast<u128>(q) * static_cast<u128>(cnt);
    }

    return static_cast<u64>(static_cast<u128>(n) - non_smooth);
}

bool is_smooth_bruteforce(u64 n) {
    if (n == 1ULL) {
        return true;
    }
    const long double root = std::sqrt(static_cast<long double>(n));
    u64 x = n;
    for (u64 p = 2; p * p <= x; ++p) {
        if (x % p != 0ULL) {
            continue;
        }
        if (static_cast<long double>(p) >= root) {
            return false;
        }
        while (x % p == 0ULL) {
            x /= p;
        }
    }
    if (x > 1ULL && static_cast<long double>(x) >= root) {
        return false;
    }
    return true;
}

u64 brute_count(u64 n) {
    u64 count = 0;
    for (u64 i = 1; i <= n; ++i) {
        if (is_smooth_bruteforce(i)) {
            ++count;
        }
    }
    return count;
}

}  // namespace

int main() {
    PrimeCounting pc;

    assert(count_smooth(100ULL, pc) == 29ULL);
    assert(count_smooth(1'000ULL, pc) == brute_count(1'000ULL));

    std::cout << count_smooth(10'000'000'000ULL, pc) << "\n";
    return 0;
}

Python

def solve_case(n):
    r = int(n**0.5)
    V = [n // i for i in range(1, r + 1)]
    for i in range(V[-1] - 1, 0, -1):
        V.append(i)
        
    S = {i: i - 1 for i in V}
    S[0] = 0
    
    for p in range(2, r + 1):
        if S[p] > S[p - 1]:
            sp = S[p - 1]
            p2 = p * p
            for v in V:
                if v < p2: break
                S[v] -= S[v // p] - sp
                
    non_smooth = 0
    for q in range(1, r + 1):
        non_smooth += S[n // q] - S[q - 1]
        
    return n - non_smooth

def solve():
    ans = solve_case(10000000000)
    return str(ans)

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

Java

public class Euler668 {

    static int getIdx(long val, long n, int r, int vCount) {
        if (val <= r) {
            return vCount - (int) val;
        }
        return (int) (n / val) - 1;
    }

    static long solveCase(long n) {
        int r = (int) Math.sqrt(n);
        long[] V = new long[2 * r + 2];
        int vCount = 0;

        for (int i = 1; i <= r; ++i) {
            V[vCount++] = n / i;
        }
        for (long i = V[vCount - 1] - 1; i > 0; --i) {
            V[vCount++] = i;
        }

        long[] S = new long[vCount];
        for (int i = 0; i < vCount; ++i) {
            S[i] = V[i] - 1;
        }

        for (long p = 2; p <= r; ++p) {
            int idxP = getIdx(p, n, r, vCount);
            int idxPminus1 = getIdx(p - 1, n, r, vCount);

            if (S[idxP] > S[idxPminus1]) {
                long sp = S[idxPminus1];
                long p2 = p * p;

                for (int i = 0; i < vCount; ++i) {
                    long v = V[i];
                    if (v < p2)
                        break;

                    long nextV = v / p;
                    int nextIdx = getIdx(nextV, n, r, vCount);
                    S[i] -= (S[nextIdx] - sp);
                }
            }
        }

        long nonSmooth = 0;
        for (int q = 1; q <= r; ++q) {
            long rightVal = n / q;
            long leftVal = q - 1;

            int rightIdx = getIdx(rightVal, n, r, vCount);
            long leftS = leftVal == 0 ? 0 : S[getIdx(leftVal, n, r, vCount)];

            nonSmooth += S[rightIdx] - leftS;
        }

        return n - nonSmooth;
    }

    public static String solve() {
        long ans = solveCase(10000000000L);
        return Long.toString(ans);
    }

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