Problem 652: Distinct Values of a Proto-logarithmic Function

Project Euler Problem 652 Solution

EulerSolve provides an optimized solution for Project Euler Problem 652, Distinct Values of a Proto-logarithmic Function, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For every ordered pair $(x,y)$ with $2 \le x,y \le N$, the problem asks for the number $D(N)$ of distinct proto-logarithmic values that can occur. The central observation used by the C++, Python, and Java implementations is that every integer $n \ge 2$ has a unique canonical form $$n=r^k,$$ where $r$ is not a perfect power and $k \ge 1$ is maximal. Once both inputs are rewritten in that form, all duplicates are controlled by primitive bases and by reduced exponent ratios, so the count can be expressed using only $L=\lfloor \log_2 N \rfloor$ and some short arithmetic tables. Mathematical Approach We write the final answer as $$D(N)=\operatorname{Rat}(N)+\operatorname{Irr}(N),$$ where the first term counts values from the equal-base branch and the second counts values from the distinct-base branch. Step 1: Rewrite Every Integer in Canonical Perfect-Power Form Every integer $n \ge 2$ can be written uniquely as $n=r^k$ with $r$ not a perfect power. Call such an $r$ a primitive base ....

Detailed mathematical approach

Problem Summary

For every ordered pair $(x,y)$ with $2 \le x,y \le N$, the problem asks for the number $D(N)$ of distinct proto-logarithmic values that can occur. The central observation used by the C++, Python, and Java implementations is that every integer $n \ge 2$ has a unique canonical form

$$n=r^k,$$

where $r$ is not a perfect power and $k \ge 1$ is maximal. Once both inputs are rewritten in that form, all duplicates are controlled by primitive bases and by reduced exponent ratios, so the count can be expressed using only $L=\lfloor \log_2 N \rfloor$ and some short arithmetic tables.

Mathematical Approach

We write the final answer as

$$D(N)=\operatorname{Rat}(N)+\operatorname{Irr}(N),$$

where the first term counts values from the equal-base branch and the second counts values from the distinct-base branch.

Step 1: Rewrite Every Integer in Canonical Perfect-Power Form

Every integer $n \ge 2$ can be written uniquely as $n=r^k$ with $r$ not a perfect power. Call such an $r$ a primitive base. For a fixed primitive base $r$, define its exponent capacity by

$$A(r)=\max\{e \ge 1 : r^e \le N\}.$$

So the numbers generated by $r$ inside the range $[2,N]$ are exactly

$$r,r^2,\dots,r^{A(r)}.$$

If we rewrite an ordered pair as

$$x=r^u,\qquad y=s^v,$$

with $1 \le u \le A(r)$ and $1 \le v \le A(s)$, then the problem depends only on the ordered primitive bases $(r,s)$ and on the reduced fraction $v/u$. This mirrors the usual identity

$$\log_{r^u}(s^v)=\frac{v}{u}\cdot\frac{\log s}{\log r},$$

which explains why the equal-base branch becomes rational and the distinct-base branch behaves like an irrational family indexed by $(r,s)$.

Step 2: Count Primitive Bases with Möbius Inversion

Let

$$U(M)=\#\{r \in \mathbb{Z}_{\ge 2} : r \le M,\ \forall a,b \in \mathbb{Z}_{\ge 2},\ r \ne a^b\}.$$

Every integer in $[2,M]$ has a unique representation $r^e$ with primitive base $r$, so counting all integers up to $M$ gives

$$M-1=\sum_{e=1}^{\lfloor \log_2 M \rfloor} U\!\left(\left\lfloor M^{1/e}\right\rfloor\right).$$

Applying Möbius inversion yields the exact formula

$$U(M)=\sum_{d=1}^{\lfloor \log_2 M \rfloor} \mu(d)\left(\left\lfloor M^{1/d}\right\rfloor-1\right).$$

Now define $B_A(N)$ to be the number of primitive bases whose exponent capacity is exactly $A$. Then

$$B_A(N)=U\!\left(\left\lfloor N^{1/A}\right\rfloor\right)-U\!\left(\left\lfloor N^{1/(A+1)}\right\rfloor\right).$$

This partitions all primitive bases according to the largest exponent they can support inside the range.

Step 3: Count Reduced Exponent Ratios

For fixed bounds $A$ and $B$, define

$$Q(A,B)=\#\{(u,v): 1 \le u \le A,\ 1 \le v \le B,\ \gcd(u,v)=1\}.$$

This is exactly the number of distinct reduced fractions $v/u$ that can appear when the denominator exponent is at most $A$ and the numerator exponent is at most $B$.

Using the standard Möbius count of visible lattice points,

$$Q(A,B)=\sum_{d=1}^{\min(A,B)} \mu(d)\left\lfloor \frac{A}{d} \right\rfloor \left\lfloor \frac{B}{d} \right\rfloor.$$

The implementations precompute this entire table once for all $1 \le A,B \le L$, where

$$L=\lfloor \log_2 N \rfloor.$$

Step 4: Rational Branch

When the primitive bases coincide, say $r=s$, the value depends only on the reduced fraction $v/u$. Therefore the total number of distinct rational values is exactly the number of reduced positive fractions with both numerator and denominator at most $L$:

$$\operatorname{Rat}(N)=Q(L,L).$$

This works because the primitive base $2$ alone already supports exponents $1,2,\dots,L$, since $2^L \le N$ by definition of $L$.

Step 5: Irrational Branch

When $r \ne s$, the value belongs to the distinct-base branch. For a fixed ordered pair of primitive bases $(r,s)$, the only collisions come from reducing the fraction $v/u$, so that pair contributes exactly $Q(A(r),A(s))$ distinct values.

How many ordered primitive-base pairs have capacities $A$ and $B$? If $A \ne B$, the capacities already force different bases, so the count is

$$P_{A,B}(N)=B_A(N)\,B_B(N).$$

If $A=B$, we must exclude choosing the same primitive base twice, so

$$P_{A,A}(N)=B_A(N)\bigl(B_A(N)-1\bigr).$$

Hence the entire irrational branch is

$$\operatorname{Irr}(N)=\sum_{A=1}^{L}\sum_{B=1}^{L} P_{A,B}(N)\,Q(A,B).$$

Worked Example: $N=5$

Here $L=\lfloor \log_2 5 \rfloor=2$. The primitive bases not exceeding $5$ are $2,3,5$. Their exponent capacities are

$$A(2)=2,\qquad A(3)=1,\qquad A(5)=1.$$

So the capacity classes are

$$B_1(5)=2,\qquad B_2(5)=1.$$

The reduced-ratio counts are

$$Q(1,1)=1,\qquad Q(1,2)=2,\qquad Q(2,1)=2,\qquad Q(2,2)=3.$$

The rational branch contributes

$$\operatorname{Rat}(5)=Q(2,2)=3,$$

corresponding to the reduced fractions $1$, $2$, and $1/2$.

For the irrational branch, the ordered-pair counts are

$$P_{1,1}(5)=2,\qquad P_{1,2}(5)=2,\qquad P_{2,1}(5)=2,\qquad P_{2,2}(5)=0.$$

Therefore

$$\operatorname{Irr}(5)=2 \cdot 1 + 2 \cdot 2 + 2 \cdot 2 = 10,$$

and the total is

$$D(5)=3+10=13.$$

This matches the first exact checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations first compute $L=\lfloor \log_2 N \rfloor$ and build the Möbius values $\mu(1),\dots,\mu(L)$ with a sieve. Because all later formulas depend only on exponent bounds up to $L$, this immediately compresses the problem from a search over numbers up to $N$ into a search over exponents up to roughly $60$ for the target input.

Next, the implementation evaluates exact integer $k$-th roots such as $\lfloor N^{1/k} \rfloor$. It starts from a floating-point estimate and then corrects it with integer power checks, so no rounding error leaks into the Möbius sums. Those exact roots feed the formula for $U(M)$, and from there the code derives every $B_A(N)$.

After that, the implementation fills the table $Q(A,B)$ for all $1 \le A,B \le L$ using the Möbius formula for coprime exponent pairs. Once this table exists, the rational part is just the single lookup $Q(L,L)$.

The remaining work is the double sum for the irrational part. Off the diagonal $A \ne B$, the number of ordered base pairs is $B_A(N)B_B(N)$. On the diagonal, the code uses $B_A(N)(B_A(N)-1)$ so that a primitive base is not paired with itself in the distinct-base branch.

For the full target $N=10^{18}$, the answer is reported modulo $10^9$. For smaller checkpoints, the implementations also support exact arithmetic, and the built-in validations verify values such as $D(5)=13$, $D(10)=69$, $D(100)=9607$, and $D(10000)=99959605$.

Complexity Analysis

Let $L=\lfloor \log_2 N \rfloor$. The Möbius sieve is $O(L)$. Building the coprime-ratio table requires

$$\sum_{A=1}^{L}\sum_{B=1}^{L} O(\min(A,B))=O(L^3)$$

operations and is the dominant asymptotic cost in the local implementations. The final irrational double sum is $O(L^2)$, and the storage cost is $O(L^2)$ because the $Q(A,B)$ table is kept in memory. For the target input $N=10^{18}$, we only have $L=59$, so the method is easily fast enough.

Footnotes and References

Problem page: Project Euler 652
Perfect powers: Wikipedia - Perfect power
Möbius function: Wikipedia - Möbius function
Möbius inversion formula: Wikipedia - Möbius inversion formula
Coprime integers: Wikipedia - Coprime integers

Problem 652 source code

C++

#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <thread>
#include <vector>
#include <functional>

using namespace std;

namespace {

constexpr uint64_t kTargetN = 1'000'000'000'000'000'000ULL;
constexpr uint64_t kMod = 1'000'000'000ULL;

int floor_log2_u64(uint64_t x) {
    int r = 0;
    while (x > 1) {
        x >>= 1;
        ++r;
    }
    return r;
}

bool pow_leq(uint64_t base, int exp, uint64_t limit) {
    unsigned __int128 res = 1;
    for (int i = 0; i < exp; ++i) {
        res *= base;
        if (res > limit) return false;
    }
    return true;
}

uint64_t int_nth_root(uint64_t n, int k) {
    if (k == 1 || n <= 1) return n;
    long double approx = std::pow(static_cast<long double>(n), 1.0L / k);
    uint64_t r = static_cast<uint64_t>(approx);
    if (r < 1) r = 1;
    while (pow_leq(r + 1, k, n)) ++r;
    while (!pow_leq(r, k, n)) --r;
    return r;
}

void mobius_sieve(int nmax, vector<int8_t>& mu) {
    mu.assign(nmax + 1, 0);
    vector<int> primes;
    primes.reserve(nmax / 2);
    mu[1] = 1;
    vector<int> spf(nmax + 1, 0);
    for (int i = 2; i <= nmax; ++i) {
        if (spf[i] == 0) {
            spf[i] = i;
            primes.push_back(i);
            mu[i] = -1;
        }
        for (int p : primes) {
            long long v = 1LL * p * i;
            if (v > nmax) break;
            spf[static_cast<int>(v)] = p;
            if (i % p == 0) {
                mu[static_cast<int>(v)] = 0;
                break;
            }
            mu[static_cast<int>(v)] = static_cast<int8_t>(-mu[i]);
        }
    }
}

// Count integers r in [2, M] that are not perfect powers.
uint64_t power_free_count_upto(uint64_t M, const vector<int8_t>& mu) {
    if (M < 2) return 0;
    int L = floor_log2_u64(M);
    int max_d = min(L, static_cast<int>(mu.size()) - 1);
    int64_t total = 0;
    for (int d = 1; d <= max_d; ++d) {
        const int8_t md = mu[d];
        if (md == 0) continue;
        uint64_t root = int_nth_root(M, d);
        if (root >= 2) {
            total += static_cast<int64_t>(md) * static_cast<int64_t>(root - 1);
        }
    }
    return static_cast<uint64_t>(total);
}

vector<vector<int64_t>> build_coprime_table(int L, const vector<int8_t>& mu) {
    vector<vector<int64_t>> F(L + 1, vector<int64_t>(L + 1, 0));
    for (int A = 1; A <= L; ++A) {
        for (int B = 1; B <= L; ++B) {
            int m = min(A, B);
            int64_t total = 0;
            for (int d = 1; d <= m; ++d) {
                int8_t md = mu[d];
                if (md == 0) continue;
                total += static_cast<int64_t>(md) * (A / d) * (B / d);
            }
            F[A][B] = total;
        }
    }
    return F;
}

vector<uint64_t> build_countA(uint64_t N, const vector<int8_t>& mu, int L) {
    vector<uint64_t> roots(L + 2, 0);
    for (int k = 1; k <= L; ++k) {
        roots[k] = int_nth_root(N, k);
    }
    roots[L + 1] = 1;

    vector<uint64_t> pf(L + 2, 0);
    for (int k = 1; k <= L + 1; ++k) {
        pf[k] = power_free_count_upto(roots[k], mu);
    }

    vector<uint64_t> countA(L + 1, 0);
    for (int A = 1; A <= L; ++A) {
        countA[A] = pf[A] - pf[A + 1];
    }
    return countA;
}

uint64_t mod_mul(uint64_t a, uint64_t b, uint64_t mod) {
    return static_cast<uint64_t>((static_cast<unsigned __int128>(a) * b) % mod);
}

unsigned __int128 compute_irrational_exact(const vector<uint64_t>& countA,
                                           const vector<vector<int64_t>>& coprime,
                                           unsigned threads) {
    const int L = static_cast<int>(countA.size()) - 1;
    if (threads <= 1 || L <= 1) {
        unsigned __int128 irr = 0;
        for (int A = 1; A <= L; ++A) {
            uint64_t ca = countA[A];
            if (ca == 0) continue;
            for (int B = 1; B <= L; ++B) {
                uint64_t cb = countA[B];
                if (cb == 0) continue;
                unsigned __int128 pairs = 0;
                if (A == B) {
                    if (ca < 2) continue;
                    pairs = static_cast<unsigned __int128>(ca) * (ca - 1);
                } else {
                    pairs = static_cast<unsigned __int128>(ca) * cb;
                }
                irr += pairs * static_cast<unsigned __int128>(coprime[A][B]);
            }
        }
        return irr;
    }

    threads = min<unsigned>(threads, static_cast<unsigned>(L));
    vector<unsigned __int128> partial(threads, 0);
    atomic<int> nextA(1);
    const int chunk = 2;

    auto worker = [&](unsigned idx) {
        unsigned __int128 local = 0;
        while (true) {
            int start = nextA.fetch_add(chunk);
            if (start > L) break;
            int end = min(L, start + chunk - 1);
            for (int A = start; A <= end; ++A) {
                uint64_t ca = countA[A];
                if (ca == 0) continue;
                for (int B = 1; B <= L; ++B) {
                    uint64_t cb = countA[B];
                    if (cb == 0) continue;
                    unsigned __int128 pairs = 0;
                    if (A == B) {
                        if (ca < 2) continue;
                        pairs = static_cast<unsigned __int128>(ca) * (ca - 1);
                    } else {
                        pairs = static_cast<unsigned __int128>(ca) * cb;
                    }
                    local += pairs * static_cast<unsigned __int128>(coprime[A][B]);
                }
            }
        }
        partial[idx] = local;
    };

    vector<thread> pool;
    pool.reserve(threads);
    for (unsigned t = 0; t < threads; ++t) {
        pool.emplace_back(worker, t);
    }
    for (auto& th : pool) th.join();

    unsigned __int128 total = 0;
    for (auto v : partial) total += v;
    return total;
}

uint64_t compute_irrational_mod(const vector<uint64_t>& countA,
                                const vector<vector<int64_t>>& coprime,
                                unsigned threads,
                                uint64_t mod) {
    const int L = static_cast<int>(countA.size()) - 1;
    if (threads == 0) threads = 1;
    if (threads <= 1 || L <= 1) {
        uint64_t irr = 0;
        for (int A = 1; A <= L; ++A) {
            uint64_t ca = countA[A];
            if (ca == 0) continue;
            uint64_t ca_mod = ca % mod;
            for (int B = 1; B <= L; ++B) {
                uint64_t cb = countA[B];
                if (cb == 0) continue;
                uint64_t pairs = 0;
                if (A == B) {
                    if (ca < 2) continue;
                    pairs = mod_mul(ca_mod, (ca - 1) % mod, mod);
                } else {
                    pairs = mod_mul(ca_mod, cb % mod, mod);
                }
                uint64_t term = mod_mul(pairs,
                                        static_cast<uint64_t>(coprime[A][B]) % mod,
                                        mod);
                irr += term;
                if (irr >= mod) irr %= mod;
            }
        }
        return irr % mod;
    }

    threads = min<unsigned>(threads, static_cast<unsigned>(L));
    vector<uint64_t> partial(threads, 0);
    atomic<int> nextA(1);
    const int chunk = 2;

    auto worker = [&](unsigned idx) {
        uint64_t local = 0;
        while (true) {
            int start = nextA.fetch_add(chunk);
            if (start > L) break;
            int end = min(L, start + chunk - 1);
            for (int A = start; A <= end; ++A) {
                uint64_t ca = countA[A];
                if (ca == 0) continue;
                uint64_t ca_mod = ca % mod;
                for (int B = 1; B <= L; ++B) {
                    uint64_t cb = countA[B];
                    if (cb == 0) continue;
                    uint64_t pairs = 0;
                    if (A == B) {
                        if (ca < 2) continue;
                        pairs = mod_mul(ca_mod, (ca - 1) % mod, mod);
                    } else {
                        pairs = mod_mul(ca_mod, cb % mod, mod);
                    }
                    uint64_t term = mod_mul(pairs,
                                            static_cast<uint64_t>(coprime[A][B]) % mod,
                                            mod);
                    local += term;
                    if (local >= mod) local %= mod;
                }
            }
        }
        partial[idx] = local % mod;
    };

    vector<thread> pool;
    pool.reserve(threads);
    for (unsigned t = 0; t < threads; ++t) {
        pool.emplace_back(worker, t);
    }
    for (auto& th : pool) th.join();

    uint64_t total = 0;
    for (uint64_t v : partial) {
        total += v;
        if (total >= mod) total %= mod;
    }
    return total % mod;
}

unsigned __int128 compute_D_exact(uint64_t N, unsigned threads) {
    const int L = floor_log2_u64(N);
    vector<int8_t> mu;
    mobius_sieve(L, mu);
    auto coprime = build_coprime_table(L, mu);
    auto countA = build_countA(N, mu, L);

    unsigned __int128 irrational = compute_irrational_exact(countA, coprime, threads);
    unsigned __int128 rational = static_cast<unsigned __int128>(coprime[L][L]);
    return rational + irrational;
}

uint64_t compute_D_mod(uint64_t N, unsigned threads, uint64_t mod) {
    const int L = floor_log2_u64(N);
    vector<int8_t> mu;
    mobius_sieve(L, mu);
    auto coprime = build_coprime_table(L, mu);
    auto countA = build_countA(N, mu, L);

    uint64_t irrational = compute_irrational_mod(countA, coprime, threads, mod);
    uint64_t rational = static_cast<uint64_t>(coprime[L][L]) % mod;
    return (rational + irrational) % mod;
}

string to_string_u128(unsigned __int128 x) {
    if (x == 0) return "0";
    string s;
    while (x > 0) {
        int digit = static_cast<int>(x % 10);
        s.push_back(static_cast<char>('0' + digit));
        x /= 10;
    }
    reverse(s.begin(), s.end());
    return s;
}

bool run_validations(unsigned threads) {
    struct TestCase {
        uint64_t n;
        uint64_t expected;
    };
    const TestCase tests[] = {
        {5, 13},
        {10, 69},
        {100, 9607},
        {10000, 99959605},
    };

    for (const auto& tc : tests) {
        unsigned __int128 value = compute_D_exact(tc.n, threads);
        uint64_t got = static_cast<uint64_t>(value);
        if (got != tc.expected) {
            cerr << "Validation failed: D(" << tc.n << ") = " << got
                 << " (expected " << tc.expected << ")\n";
            return false;
        }
    }
    cerr << "Validation checkpoints passed.\n";
    return true;
}

} // namespace

int main(int argc, char** argv) {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    uint64_t N = kTargetN;
    unsigned threads = thread::hardware_concurrency();
    if (threads == 0) threads = 1;
    bool validate = true;

    if (argc >= 2) N = stoull(argv[1]);
    if (argc >= 3) threads = max(1u, static_cast<unsigned>(stoul(argv[2])));
    if (argc >= 4) validate = (stoi(argv[3]) != 0);

    if (validate && !run_validations(min(threads, 4u))) {
        return 1;
    }

    if (N == kTargetN) {
        uint64_t answer = compute_D_mod(N, threads, kMod);
        cout << setfill('0') << setw(9) << answer << "\n";
    } else {
        unsigned __int128 exact = compute_D_exact(N, threads);
        cout << to_string_u128(exact) << "\n";
    }
    return 0;
}

Python

import math
import sys

def floor_log2(x):
    r = 0
    while x > 1:
        x >>= 1
        r += 1
    return r

def pow_leq(base, exp, limit):
    res = 1
    for _ in range(exp):
        res *= base
        if res > limit: return False
    return True

def int_nth_root(n, k):
    if k == 1 or n <= 1: return n
    approx = int(math.pow(n, 1.0 / k))
    r = max(1, approx)
    while pow_leq(r + 1, k, n): r += 1
    while not pow_leq(r, k, n): r -= 1
    return r

def mobius_sieve(nmax):
    mu = [0] * (nmax + 1)
    primes = []
    mu[1] = 1
    spf = [0] * (nmax + 1)
    for i in range(2, nmax + 1):
        if spf[i] == 0:
            spf[i] = i
            primes.append(i)
            mu[i] = -1
        for p in primes:
            v = p * i
            if v > nmax: break
            spf[v] = p
            if i % p == 0:
                mu[v] = 0
                break
            mu[v] = -mu[i]
    return mu

def power_free_count_upto(M, mu):
    if M < 2: return 0
    L = floor_log2(M)
    max_d = min(L, len(mu) - 1)
    total = 0
    for d in range(1, max_d + 1):
        md = mu[d]
        if md == 0: continue
        root = int_nth_root(M, d)
        if root >= 2:
            total += md * (root - 1)
    return total

def build_coprime_table(L, mu):
    F = [[0] * (L + 1) for _ in range(L + 1)]
    for A in range(1, L + 1):
        for B in range(1, L + 1):
            m = min(A, B)
            total = 0
            for d in range(1, m + 1):
                md = mu[d]
                if md == 0: continue
                total += md * (A // d) * (B // d)
            F[A][B] = total
    return F

def build_countA(N, mu, L):
    roots = [0] * (L + 2)
    for k in range(1, L + 1):
        roots[k] = int_nth_root(N, k)
    roots[L + 1] = 1

    pf = [0] * (L + 2)
    for k in range(1, L + 2):
        pf[k] = power_free_count_upto(roots[k], mu)

    countA = [0] * (L + 1)
    for A in range(1, L + 1):
        countA[A] = pf[A] - pf[A + 1]
    return countA

def compute_D_mod(N, mod):
    L = floor_log2(N)
    mu = mobius_sieve(L)
    coprime = build_coprime_table(L, mu)
    countA = build_countA(N, mu, L)

    irrational = 0
    for A in range(1, L + 1):
        ca = countA[A]
        if ca == 0: continue
        ca_mod = ca % mod
        for B in range(1, L + 1):
            cb = countA[B]
            if cb == 0: continue
            if A == B:
                if ca < 2: continue
                pairs = (ca_mod * ((ca - 1) % mod)) % mod
            else:
                pairs = (ca_mod * (cb % mod)) % mod
            term = (pairs * (coprime[A][B] % mod)) % mod
            irrational = (irrational + term) % mod
            
    rational = coprime[L][L] % mod
    return (rational + irrational) % mod

def solve():
    N = 1000000000000000000
    mod = 1000000000
    ans = compute_D_mod(N, mod)
    return "{:09d}".format(ans)

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

Java

public class Euler652 {
    static final long kTargetN = 1000000000000000000L;
    static final long kMod = 1000000000L;

    static int floorLog2(long x) {
        int r = 0;
        while (x > 1) {
            x >>= 1;
            ++r;
        }
        return r;
    }

    static boolean powLeq(long base, int exp, long limit) {
        long res = 1;
        // avoid overflow carefully
        for (int i = 0; i < exp; ++i) {
            // maxLimit / base
            if (res > (limit / base))
                return false;
            long test = res * base;
            if (test > limit || test < 0)
                return false;
            res = test;
        }
        return true;
    }

    static long intNthRoot(long n, int k) {
        if (k == 1 || n <= 1)
            return n;
        long approx = (long) Math.pow((double) n, 1.0 / k);
        long r = Math.max(1L, approx);
        while (powLeq(r + 1, k, n))
            ++r;
        while (!powLeq(r, k, n))
            --r;
        return r;
    }

    static int[] mobiusSieve(int nmax) {
        int[] mu = new int[nmax + 1];
        java.util.ArrayList<Integer> primes = new java.util.ArrayList<>();
        mu[1] = 1;
        int[] spf = new int[nmax + 1];
        for (int i = 2; i <= nmax; ++i) {
            if (spf[i] == 0) {
                spf[i] = i;
                primes.add(i);
                mu[i] = -1;
            }
            for (int p : primes) {
                long v = (long) p * i;
                if (v > nmax)
                    break;
                spf[(int) v] = p;
                if (i % p == 0) {
                    mu[(int) v] = 0;
                    break;
                }
                mu[(int) v] = -mu[i];
            }
        }
        return mu;
    }

    static long powerFreeCountUpto(long M, int[] mu) {
        if (M < 2)
            return 0;
        int L = floorLog2(M);
        int max_d = Math.min(L, mu.length - 1);
        long total = 0;
        for (int d = 1; d <= max_d; ++d) {
            int md = mu[d];
            if (md == 0)
                continue;
            long root = intNthRoot(M, d);
            if (root >= 2) {
                total += (long) md * (root - 1);
            }
        }
        return total;
    }

    static long[][] buildCoprimeTable(int L, int[] mu) {
        long[][] F = new long[L + 1][L + 1];
        for (int A = 1; A <= L; ++A) {
            for (int B = 1; B <= L; ++B) {
                int m = Math.min(A, B);
                long total = 0;
                for (int d = 1; d <= m; ++d) {
                    int md = mu[d];
                    if (md == 0)
                        continue;
                    total += (long) md * (A / d) * (B / d);
                }
                F[A][B] = total;
            }
        }
        return F;
    }

    static long[] buildCountA(long N, int[] mu, int L) {
        long[] roots = new long[L + 2];
        for (int k = 1; k <= L; ++k) {
            roots[k] = intNthRoot(N, k);
        }
        roots[L + 1] = 1;

        long[] pf = new long[L + 2];
        for (int k = 1; k <= L + 1; ++k) {
            pf[k] = powerFreeCountUpto(roots[k], mu);
        }

        long[] countA = new long[L + 1];
        for (int A = 1; A <= L; ++A) {
            countA[A] = pf[A] - pf[A + 1];
        }
        return countA;
    }

    static long computeDMod(long N, long mod) {
        int L = floorLog2(N);
        int[] mu = mobiusSieve(L);
        long[][] coprime = buildCoprimeTable(L, mu);
        long[] countA = buildCountA(N, mu, L);

        long irrational = 0;
        for (int A = 1; A <= L; ++A) {
            long ca = countA[A];
            if (ca == 0)
                continue;
            long ca_mod = ca % mod;
            for (int B = 1; B <= L; ++B) {
                long cb = countA[B];
                if (cb == 0)
                    continue;
                long pairs = 0;
                if (A == B) {
                    if (ca < 2)
                        continue;
                    pairs = (ca_mod * ((ca - 1) % mod)) % mod;
                } else {
                    pairs = (ca_mod * (cb % mod)) % mod;
                }
                long term = (pairs * (coprime[A][B] % mod)) % mod;
                irrational = (irrational + term) % mod;
            }
        }
        long rational = coprime[L][L] % mod;
        return (rational + irrational) % mod;
    }

    public static String solve() {
        long N = kTargetN;
        long ans = computeDMod(N, kMod);
        return String.format("%09d", ans);
    }

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