Problem 639: Summing a Multiplicative Function

Project Euler Problem 639 Solution

EulerSolve provides an optimized solution for Project Euler Problem 639, Summing a Multiplicative Function, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each positive integer $k$, define $$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k,$$ where $\operatorname{rad}(n)$ is the product of the distinct prime divisors of $n$. The problem asks for $$\sum_{k=1}^{50} S_k(10^{12}) \pmod{10^9+7}.$$ The direct approach would require evaluating $\operatorname{rad}(n)$ for every $n\le 10^{12}$, which is hopeless. The implementation instead starts from the easier power sum $\sum n^k$, then corrects it prime by prime until the local behavior matches $\operatorname{rad}(n)^k$. Mathematical Approach Fix one value of $k$. The task is to compute $S_k(N)$, and then repeat that computation for $k=1,2,\dots,50$. Step 1: Start from the easier sum $\sum n^k$ The multiplicative function $n^k$ has local values $$n^k=\prod_{p^a\parallel n} p^{ak}.$$ By contrast, $$\operatorname{rad}(n)^k=\prod_{p^a\parallel n} p^k.$$ So the only difference is what happens when a prime appears with exponent $a\ge 2$. If a prime occurs only once, both functions contribute the same factor $p^k$. The implementation therefore begins with the summatory function $$P_k(x)=\sum_{n\le x} n^k,$$ which can be evaluated quickly for many different $x$ by a Faulhaber polynomial modulo $10^9+7$....

Detailed mathematical approach

Problem Summary

For each positive integer $k$, define

$$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k,$$

where $\operatorname{rad}(n)$ is the product of the distinct prime divisors of $n$. The problem asks for

$$\sum_{k=1}^{50} S_k(10^{12}) \pmod{10^9+7}.$$

The direct approach would require evaluating $\operatorname{rad}(n)$ for every $n\le 10^{12}$, which is hopeless. The implementation instead starts from the easier power sum $\sum n^k$, then corrects it prime by prime until the local behavior matches $\operatorname{rad}(n)^k$.

Mathematical Approach

Fix one value of $k$. The task is to compute $S_k(N)$, and then repeat that computation for $k=1,2,\dots,50$.

Step 1: Start from the easier sum $\sum n^k$

The multiplicative function $n^k$ has local values

$$n^k=\prod_{p^a\parallel n} p^{ak}.$$

By contrast,

$$\operatorname{rad}(n)^k=\prod_{p^a\parallel n} p^k.$$

So the only difference is what happens when a prime appears with exponent $a\ge 2$. If a prime occurs only once, both functions contribute the same factor $p^k$.

The implementation therefore begins with the summatory function

$$P_k(x)=\sum_{n\le x} n^k,$$

which can be evaluated quickly for many different $x$ by a Faulhaber polynomial modulo $10^9+7$.

Step 2: Correct one prime at a time

Suppose we have a temporary multiplicative weight in which the primes already processed behave like $\operatorname{rad}(n)^k$, while the remaining primes still behave like $n^k$. For the next prime $p$, only prime powers $p^a$ with $a\ge 2$ need correction.

For one local factor we want to replace

$$p^{ak}\quad\text{by}\quad p^k \qquad (a\ge 2).$$

The difference is

$$p^{ak}-p^k=p^k\left(p^{(a-1)k}-1\right).$$

Using a geometric series, this becomes

$$p^{ak}-p^k=p^k(p^k-1)\left(1+p^k+p^{2k}+\cdots+p^{(a-2)k}\right).$$

Therefore the update coefficient for prime $p$ is

$$t_p=p^k(p^k-1).$$

Step 3: Derive the summatory recurrence

Let $F(x)$ be the current summatory function before correcting prime $p$. Then subtracting

$$t_p\sum_{a\ge 2} F\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right)$$

exactly changes every local factor $p^{ak}$ with $a\ge 2$ into $p^k$.

To see why, consider an integer $n=p^e m$ with $p\nmid m$. Its old local contribution at $p$ is $p^{ek}$. The subtraction contributes

$$t_p\sum_{a=2}^{e} p^{(e-a)k}=p^k(p^k-1)\sum_{j=0}^{e-2} p^{jk}=p^{ek}-p^k.$$

After subtraction, the new local contribution is therefore

$$p^{ek}-(p^{ek}-p^k)=p^k,$$

which is exactly what $\operatorname{rad}(p^e)^k$ requires. If $e=0$ or $e=1$, no subtraction occurs, and the factor is already correct.

So the prime-by-prime transition is

$$F_{\text{new}}(x)=F_{\text{old}}(x)-t_p\sum_{a\ge 2} F_{\text{old}}\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right).$$

Step 4: Only primes up to $\sqrt{N}$ matter

If $p>\sqrt{N}$, then $p^2>N$. Such a prime can appear in any $n\le N$ only with exponent $0$ or $1$. But for exponents $0$ and $1$, the weights $n^k$ and $\operatorname{rad}(n)^k$ already agree. Hence corrections are needed only for primes

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

This is why the implementation sieves primes only up to $\lfloor\sqrt{N}\rfloor$.

Step 5: Compress all reachable arguments

The recurrence never asks for arbitrary values of $x$. Starting from $N$, every transition divides by some $p^a$ with $a\ge 2$. After several corrections, every queried argument has the form

$$\left\lfloor\frac{N}{m}\right\rfloor,$$

where $m$ is a powerful number, meaning that every prime exponent in $m$ is either $0$ or at least $2$.

This is the crucial compression step. Instead of storing values for all $1\le x\le N$, the implementation stores values only for the distinct quotients produced by powerful denominators. These quotients are sorted once, indexed once, and then reused for all fifty values of $k$.

Step 6: Sum over $k=1,2,\dots,50$

For each fixed $k$, the algorithm starts from $P_k(x)$, applies the prime corrections in increasing prime order, and obtains

$$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k \pmod{10^9+7}.$$

The final answer is simply

$$\sum_{k=1}^{50} S_k(N)\pmod{10^9+7}.$$

Worked Example: $S_1(10)=41$

For $k=1$, the base sum is

$$P_1(10)=1+2+\cdots+10=55.$$

Only primes up to $\sqrt{10}$ need correction, so only $p=2$ and $p=3$ matter.

For $p=2$,

$$t_2=2(2-1)=2.$$

The correction uses $2^2$ and $2^3$:

$$2\left(P_1\!\left(\left\lfloor\frac{10}{4}\right\rfloor\right)+P_1\!\left(\left\lfloor\frac{10}{8}\right\rfloor\right)\right)=2\left(P_1(2)+P_1(1)\right)=2(3+1)=8.$$

So the total becomes $55-8=47$.

For $p=3$,

$$t_3=3(3-1)=6,$$

and only $3^2$ contributes:

$$6\,P_1\!\left(\left\lfloor\frac{10}{9}\right\rfloor\right)=6\,P_1(1)=6.$$

Therefore

$$S_1(10)=55-8-6=41.$$

This matches the direct check

$$\operatorname{rad}(1),\dots,\operatorname{rad}(10)=1,2,3,2,5,6,7,2,3,10,$$

whose sum is $41$.

How the Code Works

The C++, Python, and Java implementations follow the same structure. First they generate every prime up to $\lfloor\sqrt{N}\rfloor$. Then they precompute Faulhaber-style polynomial coefficients modulo $10^9+7$, so that $\sum_{n\le x} n^k$ can be evaluated quickly for any needed quotient $x$.

Next they enumerate the powerful denominators that can arise in the recurrence. From those denominators they build the distinct quotient states $\left\lfloor N/m\right\rfloor$, sort them in descending order, and create index tables so every future division lands on an already known state.

For each $k\in\{1,\dots,50\}$, every stored state is initialized with the corresponding power sum $P_k(x)$. The implementation then processes primes in increasing order. For a given prime $p$, it computes the factor $t_p=p^k(p^k-1)$, and for every state $x\ge p^2$ it subtracts the indexed states at $\left\lfloor x/p^2\right\rfloor,\left\lfloor x/p^3\right\rfloor,\dots$.

The states are updated from larger quotients to smaller quotients. That ordering is important: when the implementation needs the value at $\left\lfloor x/p^a\right\rfloor$, it still reads the pre-update value for the current prime, which is exactly what the recurrence requires.

After all relevant primes have been processed, the state corresponding to $x=N$ is $S_k(N)$. The implementation adds that value to the running total and continues with the next $k$.

Complexity Analysis

Sieving primes up to $\sqrt{N}$ costs $O(\sqrt{N}\log\log N)$ time and $O(\sqrt{N})$ memory. The Faulhaber preprocessing for $k\le 50$ is tiny compared with the rest of the computation.

The main optimization is that the recurrence is evaluated only on compressed quotient states of the form $\left\lfloor N/m\right\rfloor$ with powerful $m$. That state set is far smaller than $\{1,2,\dots,N\}$, which is what makes $N=10^{12}$ feasible.

For each $k$, the runtime is driven by the number of reachable pairs consisting of a quotient state and a prime power $p^a$ with $a\ge 2$. In other words, the algorithm scales with the compressed transition graph rather than with all integers up to $N$. Memory is proportional to the number of stored quotient states together with the prime list and a few auxiliary tables.

Footnotes and References

Problem page: https://projecteuler.net/problem=639
Radical of an integer: Wikipedia - Radical of an integer
Faulhaber's formula: Wikipedia - Faulhaber's formula
Bernoulli numbers: Wikipedia - Bernoulli number
Powerful numbers: Wikipedia - Powerful number
Multiplicative functions: Wikipedia - Multiplicative function

Problem 639 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

constexpr int MOD = 1'000'000'007;
constexpr u64 MAIN_N = 1'000'000'000'000ULL;
constexpr int K_MAX = 50;

int mod_pow(long long a, long long e) {
    long long r = 1 % MOD;
    a %= MOD;
    while (e > 0) {
        if (e & 1LL) r = (r * a) % MOD;
        a = (a * a) % MOD;
        e >>= 1LL;
    }
    return static_cast<int>(r);
}

int add_mod(int a, int b) {
    int s = a + b;
    if (s >= MOD) s -= MOD;
    return s;
}

int sub_mod(int a, int b) {
    int s = a - b;
    if (s < 0) s += MOD;
    return s;
}

int mul_mod(long long a, long long b) {
    return static_cast<int>((static_cast<u128>(a) * b) % MOD);
}

std::vector<int> primes_up_to(int n) {
    if (n < 2) return {};
    std::vector<std::uint8_t> is_prime(static_cast<std::size_t>(n + 1), 1U);
    is_prime[0] = is_prime[1] = 0U;
    for (int i = 2; 1LL * i * i <= n; ++i) {
        if (!is_prime[static_cast<std::size_t>(i)]) continue;
        for (int j = i * i; j <= n; j += i) {
            is_prime[static_cast<std::size_t>(j)] = 0U;
        }
    }
    std::vector<int> primes;
    primes.reserve(static_cast<std::size_t>(n / 10));
    for (int i = 2; i <= n; ++i) {
        if (is_prime[static_cast<std::size_t>(i)]) primes.push_back(i);
    }
    return primes;
}

struct FaulhaberData {
    int K = 0;
    std::vector<std::vector<int>> coeff;

    explicit FaulhaberData(int max_k) : K(max_k), coeff(static_cast<std::size_t>(max_k + 1)) {
        std::vector<int> fact(static_cast<std::size_t>(K + 2), 1);
        std::vector<int> inv_fact(static_cast<std::size_t>(K + 2), 1);
        std::vector<int> inv(static_cast<std::size_t>(K + 3), 0);

        for (int i = 1; i <= K + 1; ++i) {
            fact[static_cast<std::size_t>(i)] = mul_mod(fact[static_cast<std::size_t>(i - 1)], i);
        }
        inv_fact[static_cast<std::size_t>(K + 1)] = mod_pow(fact[static_cast<std::size_t>(K + 1)], MOD - 2);
        for (int i = K + 1; i >= 1; --i) {
            inv_fact[static_cast<std::size_t>(i - 1)] =
                mul_mod(inv_fact[static_cast<std::size_t>(i)], i);
        }
        for (int i = 1; i <= K + 2; ++i) {
            inv[static_cast<std::size_t>(i)] = mod_pow(i, MOD - 2);
        }

        auto ncr = [&](int n, int r) -> int {
            if (r < 0 || r > n) return 0;
            return mul_mod(mul_mod(fact[static_cast<std::size_t>(n)], inv_fact[static_cast<std::size_t>(r)]),
                           inv_fact[static_cast<std::size_t>(n - r)]);
        };

        std::vector<int> B(static_cast<std::size_t>(K + 1), 0);
        B[0] = 1;
        if (K >= 1) B[1] = (MOD + 1) / 2;

        for (int k = 2; k <= K; k += 2) {
            int s = 0;
            for (int i = 0; i < k; ++i) {
                int t = mul_mod(ncr(k, i), B[static_cast<std::size_t>(i)]);
                t = mul_mod(t, inv[static_cast<std::size_t>(k - i + 1)]);
                s = add_mod(s, t);
            }
            B[static_cast<std::size_t>(k)] = sub_mod(1, s);
        }

        for (int k = 1; k <= K; ++k) {
            std::vector<int> v(static_cast<std::size_t>(k + 2), 0);
            for (int i = 1; i < k; ++i) {
                int t = B[static_cast<std::size_t>(k + 1 - i)];
                t = mul_mod(t, ncr(k, i));
                t = mul_mod(t, inv[static_cast<std::size_t>(k + 1 - i)]);
                v[static_cast<std::size_t>(i)] = t;
            }
            v[static_cast<std::size_t>(k)] = (MOD + 1) / 2;
            v[static_cast<std::size_t>(k + 1)] = inv[static_cast<std::size_t>(k + 1)];
            coeff[static_cast<std::size_t>(k)] = std::move(v);
        }
    }

    int eval(u64 n, int k) const {
        const auto& v = coeff[static_cast<std::size_t>(k)];
        int n_mod = static_cast<int>(n % static_cast<u64>(MOD));
        int m = 1;
        int acc = 0;
        for (int c : v) {
            acc = add_mod(acc, mul_mod(m, c));
            m = mul_mod(m, n_mod);
        }
        return acc;
    }
};

struct ValData {
    std::vector<int> primes;
    std::vector<std::vector<u64>> V;
};

void collect_powerful(std::vector<u64>& out, u64 m, int i, u64 n, const std::vector<int>& primes) {
    out.push_back(m);

    if (static_cast<u128>(primes[static_cast<std::size_t>(i)]) * m > n) {
        return;
    }

    collect_powerful(out, m * static_cast<u64>(primes[static_cast<std::size_t>(i)]), i, n, primes);

    const int lP = static_cast<int>(primes.size());
    for (int j = i + 1; j < lP; ++j) {
        const u64 p = static_cast<u64>(primes[static_cast<std::size_t>(j)]);
        const u128 mm = static_cast<u128>(m) * p * p;
        if (mm > n) {
            return;
        }
        collect_powerful(out, static_cast<u64>(mm), j, n, primes);
    }
}

ValData build_val_data(u64 n) {
    const int lim = static_cast<int>(std::sqrt(static_cast<long double>(n)));
    std::vector<int> P = primes_up_to(lim);
    const int lP = static_cast<int>(P.size());

    std::vector<std::vector<u64>> V(static_cast<std::size_t>(lP + 1));
    V[static_cast<std::size_t>(lP)] = {n};

    std::unordered_set<u64> S;
    S.reserve(1 << 16);

    for (int i = lP - 1; i >= 0; --i) {
        std::vector<u64> arr;
        arr.reserve(256);
        arr.push_back(1ULL);

        const u64 p = static_cast<u64>(P[static_cast<std::size_t>(i)]);
        collect_powerful(arr, p * p, i, n, P);

        for (u64 x : arr) {
            S.insert(n / x);
        }

        std::vector<u64> cur;
        cur.reserve(S.size());
        for (u64 x : S) {
            cur.push_back(x);
        }
        std::sort(cur.begin(), cur.end(), std::greater<u64>());
        V[static_cast<std::size_t>(i)] = std::move(cur);
    }

    return {std::move(P), std::move(V)};
}

int solve_total(u64 n, int K) {
    const FaulhaberData F(K);
    const ValData data = build_val_data(n);

    const auto& P = data.primes;
    const auto& V = data.V;
    const int lP = static_cast<int>(P.size());

    const auto& V0 = V[0];
    std::unordered_map<u64, int> idx;
    idx.reserve(V0.size() * 2 + 16);
    for (int i = 0; i < static_cast<int>(V0.size()); ++i) {
        idx.emplace(V0[static_cast<std::size_t>(i)], i);
    }

    std::vector<std::vector<int>> Vidx(static_cast<std::size_t>(lP + 1));
    for (int i = 0; i <= lP; ++i) {
        const auto& vec = V[static_cast<std::size_t>(i)];
        auto& idv = Vidx[static_cast<std::size_t>(i)];
        idv.reserve(vec.size());
        for (u64 x : vec) {
            idv.push_back(idx[x]);
        }
    }

    int total = 0;
    for (int k = 1; k <= K; ++k) {
        std::vector<int> Svals(V0.size(), 0);
        for (int id : Vidx[0]) {
            Svals[static_cast<std::size_t>(id)] = F.eval(V0[static_cast<std::size_t>(id)], k);
        }

        for (int i = 0; i < lP; ++i) {
            const u64 p = static_cast<u64>(P[static_cast<std::size_t>(i)]);
            const u64 p2 = p * p;
            const int pk = mod_pow(static_cast<int>(p % MOD), k);
            const int t = mul_mod(pk, sub_mod(pk, 1));

            const auto& vec = V[static_cast<std::size_t>(i + 1)];
            const auto& idv = Vidx[static_cast<std::size_t>(i + 1)];

            for (std::size_t pos = 0; pos < vec.size(); ++pos) {
                const u64 x = vec[pos];
                if (x < p2) break;

                const int idx_x = idv[pos];
                int cur = Svals[static_cast<std::size_t>(idx_x)];

                u64 pp = p2;
                while (pp <= x) {
                    const auto it = idx.find(x / pp);
                    cur = sub_mod(cur, mul_mod(t, Svals[static_cast<std::size_t>(it->second)]));
                    if (pp > x / p) break;
                    pp *= p;
                }
                Svals[static_cast<std::size_t>(idx_x)] = cur;
            }
        }

        total = add_mod(total, Svals[static_cast<std::size_t>(idx[n])]);
    }

    return total;
}

bool run_validations() {
    const int s1_10 = solve_total(10ULL, 1);
    if (s1_10 != 41) {
        std::cerr << "Validation failed: S_1(10)=" << s1_10 << "\n";
        return false;
    }

    const int s1_100 = solve_total(100ULL, 1);
    if (s1_100 != 3512) {
        std::cerr << "Validation failed: S_1(100)=" << s1_100 << "\n";
        return false;
    }

    const int sum2_100 = solve_total(100ULL, 2);
    const int s2_100 = sub_mod(sum2_100, s1_100);
    if (s2_100 != 208090) {
        std::cerr << "Validation failed: S_2(100)=" << s2_100 << "\n";
        return false;
    }

    const int s1_10000 = solve_total(10000ULL, 1);
    if (s1_10000 != 35252550) {
        std::cerr << "Validation failed: S_1(10000)=" << s1_10000 << "\n";
        return false;
    }

    const int sum3_1e8 = solve_total(100000000ULL, 3);
    if (sum3_1e8 != 338787512) {
        std::cerr << "Validation failed: sum_{k<=3} S_k(1e8)=" << sum3_1e8 << "\n";
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool validate = true;
    for (int i = 1; i < argc; ++i) {
        std::string arg(argv[i]);
        if (arg == "--no-validate") {
            validate = false;
        }
    }

    if (validate && !run_validations()) return 1;
    std::cout << solve_total(MAIN_N, K_MAX) << '\n';
    return 0;
}

Python

import math

MOD = 1000000007
MAIN_N = 1000000000000
K_MAX = 50

def mod_pow(a, e):
    return pow(a, e, MOD)

def add_mod(a, b):
    return (a + b) % MOD

def sub_mod(a, b):
    return (a - b) % MOD

def mul_mod(a, b):
    return (a * b) % MOD

def primes_up_to(n):
    if n < 2: return []
    is_prime = bytearray(n + 1)
    for i in range(n + 1): is_prime[i] = 1
    is_prime[0] = is_prime[1] = 0
    primes = []
    for i in range(2, n + 1):
        if not is_prime[i]: continue
        primes.append(i)
        if i * i <= n:
            for j in range(i * i, n + 1, i):
                is_prime[j] = 0
    return primes

class FaulhaberData:
    def __init__(self, max_k):
        self.K = max_k
        self.coeff = [[] for _ in range(max_k + 1)]
        fact = [1] * (max_k + 2)
        inv_fact = [1] * (max_k + 2)
        inv = [0] * (max_k + 3)
        for i in range(1, max_k + 2): fact[i] = (fact[i - 1] * i) % MOD
        inv_fact[max_k + 1] = mod_pow(fact[max_k + 1], MOD - 2)
        for i in range(max_k + 1, 0, -1):
            inv_fact[i - 1] = (inv_fact[i] * i) % MOD
        for i in range(1, max_k + 3):
            inv[i] = mod_pow(i, MOD - 2)
            
        def ncr(n, r):
            if r < 0 or r > n: return 0
            return (fact[n] * inv_fact[r] * inv_fact[n - r]) % MOD
            
        B = [0] * (max_k + 1)
        B[0] = 1
        if max_k >= 1: B[1] = (MOD + 1) // 2
        
        for k in range(2, max_k + 1, 2):
            s = 0
            for i in range(k):
                t = (ncr(k, i) * B[i]) % MOD
                t = (t * inv[k - i + 1]) % MOD
                s = (s + t) % MOD
            B[k] = (1 - s) % MOD
            
        for k in range(1, max_k + 1):
            v = [0] * (k + 2)
            for i in range(1, k):
                t = (B[k + 1 - i] * ncr(k, i)) % MOD
                t = (t * inv[k + 1 - i]) % MOD
                v[i] = t
            v[k] = (MOD + 1) // 2
            v[k + 1] = inv[k + 1]
            self.coeff[k] = v
            
    def eval_poly(self, n, k):
        v = self.coeff[k]
        n_mod = n % MOD
        m = 1
        acc = 0
        for c in v:
            acc = (acc + m * c) % MOD
            m = (m * n_mod) % MOD
        return acc

def build_val_data(n):
    lim = int(math.sqrt(n))
    P = primes_up_to(lim)
    lP = len(P)
    V = [[] for _ in range(lP + 1)]
    V[lP] = [n]
    
    S = set()
    
    def collect_powerful(out, m, i, n_val):
        out.append(m)
        if P[i] * m > n_val: return
        collect_powerful(out, m * P[i], i, n_val)
        for j in range(i + 1, lP):
            p = P[j]
            mm = m * p * p
            if mm > n_val: return
            collect_powerful(out, mm, j, n_val)
            
    for i in range(lP - 1, -1, -1):
        arr = [1]
        p = P[i]
        collect_powerful(arr, p * p, i, n)
        for x in arr:
            S.add(n // x)
            
        cur = sorted(list(S), reverse=True)
        V[i] = cur
        
    return P, V

def solve_total(n, K_val):
    F = FaulhaberData(K_val)
    P, V = build_val_data(n)
    lP = len(P)
    V0 = V[0]
    
    idx = {}
    for i, vl in enumerate(V0):
        idx[vl] = i
        
    Vidx = [[] for _ in range(lP + 1)]
    for i in range(lP + 1):
        vec = V[i]
        Vidx[i] = [idx[x] for x in vec]
        
    total = 0
    for k in range(1, K_val + 1):
        Svals = [0] * len(V0)
        for i_v in Vidx[0]:
            Svals[i_v] = F.eval_poly(V0[i_v], k)
            
        for i in range(lP):
            p = P[i]
            p2 = p * p
            pk = mod_pow(p % MOD, k)
            t = (pk * (pk - 1)) % MOD
            
            vec = V[i + 1]
            idv = Vidx[i + 1]
            
            for pos in range(len(vec)):
                x = vec[pos]
                if x < p2: break
                
                idx_x = idv[pos]
                cur = Svals[idx_x]
                
                pp = p2
                while pp <= x:
                    it_val = idx.get(x // pp)
                    if it_val is not None:
                        cur = (cur - t * Svals[it_val]) % MOD
                    if pp > x // p: break
                    pp *= p
                Svals[idx_x] = cur
                
        total = (total + Svals[idx[n]]) % MOD
        
    return total

def solve():
    return str(solve_total(MAIN_N, K_MAX))

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

Java

import java.util.*;

public class Euler639 {
    static final int MOD = 1000000007;

    static int modPow(long a, long e) {
        long r = 1 % MOD;
        a %= MOD;
        while (e > 0) {
            if ((e & 1) != 0)
                r = (r * a) % MOD;
            a = (a * a) % MOD;
            e >>= 1;
        }
        return (int) r;
    }

    static int addMod(int a, int b) {
        int s = a + b;
        if (s >= MOD)
            s -= MOD;
        return s;
    }

    static int subMod(int a, int b) {
        int s = a - b;
        if (s < 0)
            s += MOD;
        return s;
    }

    static int mulMod(long a, long b) {
        return (int) ((a * b) % MOD);
    }

    static ArrayList<Integer> primesUpTo(int n) {
        ArrayList<Integer> primes = new ArrayList<>();
        if (n < 2)
            return primes;
        boolean[] comp = new boolean[n + 1];
        comp[0] = comp[1] = true;
        for (int i = 2; i <= n; i++) {
            if (comp[i])
                continue;
            primes.add(i);
            if ((long) i * i <= n) {
                for (int j = i * i; j <= n; j += i)
                    comp[j] = true;
            }
        }
        return primes;
    }

    static class FaulhaberData {
        int K;
        int[][] coeff;

        FaulhaberData(int maxK) {
            K = maxK;
            coeff = new int[K + 1][];
            int[] fact = new int[K + 2];
            int[] invFact = new int[K + 2];
            int[] inv = new int[K + 3];

            Arrays.fill(fact, 1);
            for (int i = 1; i <= K + 1; ++i)
                fact[i] = mulMod(fact[i - 1], i);
            invFact[K + 1] = modPow(fact[K + 1], MOD - 2);
            for (int i = K + 1; i >= 1; --i)
                invFact[i - 1] = mulMod(invFact[i], i);
            for (int i = 1; i <= K + 2; ++i)
                inv[i] = modPow(i, MOD - 2);

            int[] B = new int[K + 1];
            B[0] = 1;
            if (K >= 1)
                B[1] = (MOD + 1) / 2;

            for (int k = 2; k <= K; k += 2) {
                int s = 0;
                for (int i = 0; i < k; ++i) {
                    int t = mulMod(ncr(k, i, fact, invFact), B[i]);
                    t = mulMod(t, inv[k - i + 1]);
                    s = addMod(s, t);
                }
                B[k] = subMod(1, s);
            }

            for (int k = 1; k <= K; ++k) {
                int[] v = new int[k + 2];
                for (int i = 1; i < k; ++i) {
                    int t = B[k + 1 - i];
                    t = mulMod(t, ncr(k, i, fact, invFact));
                    t = mulMod(t, inv[k + 1 - i]);
                    v[i] = t;
                }
                v[k] = (MOD + 1) / 2;
                v[k + 1] = inv[k + 1];
                coeff[k] = v;
            }
        }

        int ncr(int n, int r, int[] fact, int[] invFact) {
            if (r < 0 || r > n)
                return 0;
            return mulMod(mulMod(fact[n], invFact[r]), invFact[n - r]);
        }

        int eval(long n, int k) {
            int[] v = coeff[k];
            int nMod = (int) (n % MOD);
            int m = 1;
            int acc = 0;
            for (int c : v) {
                acc = addMod(acc, mulMod(m, c));
                m = mulMod(m, nMod);
            }
            return acc;
        }
    }

    static class ValData {
        ArrayList<Integer> primes;
        ArrayList<Long>[] V;
    }

    static void collectPowerful(ArrayList<Long> out, long m, int i, long n, ArrayList<Integer> primes) {
        out.add(m);
        if (primes.get(i) * m > n)
            return;
        collectPowerful(out, m * primes.get(i), i, n, primes);
        int lP = primes.size();
        for (int j = i + 1; j < lP; ++j) {
            long p = primes.get(j);
            long mm = m * p * p;
            if (mm > n)
                return;
            collectPowerful(out, mm, j, n, primes);
        }
    }

    @SuppressWarnings("unchecked")
    static ValData buildValData(long n) {
        int lim = (int) Math.sqrt(n);
        ArrayList<Integer> P = primesUpTo(lim);
        int lP = P.size();
        ArrayList<Long>[] V = new ArrayList[lP + 1];
        V[lP] = new ArrayList<>();
        V[lP].add(n);

        HashSet<Long> S = new HashSet<>();
        for (int i = lP - 1; i >= 0; --i) {
            ArrayList<Long> arr = new ArrayList<>();
            arr.add(1L);
            long p = P.get(i);
            collectPowerful(arr, p * p, i, n, P);
            for (long x : arr)
                S.add(n / x);
            ArrayList<Long> cur = new ArrayList<>(S);
            Collections.sort(cur, Collections.reverseOrder());
            V[i] = cur;
        }

        ValData res = new ValData();
        res.primes = P;
        res.V = V;
        return res;
    }

    static int solveTotal(long n, int K) {
        FaulhaberData F = new FaulhaberData(K);
        ValData data = buildValData(n);
        ArrayList<Integer> P = data.primes;
        ArrayList<Long>[] V = data.V;
        int lP = P.size();
        ArrayList<Long> V0 = V[0];

        HashMap<Long, Integer> idx = new HashMap<>();
        for (int i = 0; i < V0.size(); ++i) {
            idx.put(V0.get(i), i);
        }

        int[][] Vidx = new int[lP + 1][];
        for (int i = 0; i <= lP; ++i) {
            ArrayList<Long> vec = V[i];
            int[] idv = new int[vec.size()];
            for (int pos = 0; pos < vec.size(); pos++) {
                idv[pos] = idx.get(vec.get(pos));
            }
            Vidx[i] = idv;
        }

        int total = 0;
        int[] Svals = new int[V0.size()];

        for (int k = 1; k <= K; ++k) {
            for (int id : Vidx[0]) {
                Svals[id] = F.eval(V0.get(id), k);
            }

            for (int i = 0; i < lP; ++i) {
                long p = P.get(i);
                long p2 = p * p;
                int pk = modPow(p % MOD, k);
                int t = mulMod(pk, subMod(pk, 1));

                ArrayList<Long> vec = V[i + 1];
                int[] idv = Vidx[i + 1];

                for (int pos = 0; pos < vec.size(); ++pos) {
                    long x = vec.get(pos);
                    if (x < p2)
                        break;
                    int idx_x = idv[pos];
                    int cur = Svals[idx_x];

                    long pp = p2;
                    while (pp <= x) {
                        Integer it = idx.get(x / pp);
                        if (it != null) {
                            cur = subMod(cur, mulMod(t, Svals[it]));
                        }
                        if (pp > x / p)
                            break;
                        pp *= p;
                    }
                    Svals[idx_x] = cur;
                }
            }
            total = addMod(total, Svals[idx.get(n)]);
        }
        return total;
    }

    public static String solve() {
        return Integer.toString(solveTotal(1000000000000L, 50));
    }

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