Problem 484: Arithmetic Derivative

Project Euler Problem 484 Solution

EulerSolve provides an optimized solution for Project Euler Problem 484, Arithmetic Derivative, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We must evaluate $$S(N)=\sum_{n=2}^{N}\gcd(n,n'),$$ where $n'$ is the arithmetic derivative, defined by $p'=1$ for every prime $p$ and by the product rule $(ab)'=a'b+ab'$. For $n=\prod p^{e_p}$, this implies $$n'=n\sum_{p^{e_p}\parallel n}\frac{e_p}{p}.$$ The target value uses $N=5\times 10^{15}$, so factoring every integer up to $N$ is completely infeasible. The solution therefore turns $\gcd(n,n')$ into a multiplicative quantity and sums it by grouping integers according to their repeated prime factors. Mathematical Approach Write $$g(n)=\gcd(n,n').$$ The whole problem is to compute $\sum_{n\le N} g(n)$ and then exclude $n=1$. Step 1: Find the Local Prime-Power Factor Fix a prime $p$, and write $n=p^e m$ with $p\nmid m$. Using the product rule, $$n'=(p^e m)'=e\,p^{e-1}m+p^e m'=p^{e-1}(e\,m+p\,m').$$ Because $m$ is not divisible by $p$, the $p$-adic valuation of $n'$ behaves in a simple way: $$v_p(n')=\begin{cases} e-1,& p\nmid e,\\ \ge e,& p\mid e. \end{cases}$$ Therefore the $p$-part of $\gcd(n,n')$ is $$v_p(g(n))=\min\!\bigl(e,v_p(n')\bigr)=\begin{cases} e-1,& p\nmid e,\\ e,& p\mid e. \end{cases}$$ Equivalently, for one prime power, $$g(p^e)=\begin{cases} p^{e-1},& p\nmid e,\\ p^e,& p\mid e....

Detailed mathematical approach

Problem Summary

We must evaluate

$$S(N)=\sum_{n=2}^{N}\gcd(n,n'),$$

where $n'$ is the arithmetic derivative, defined by $p'=1$ for every prime $p$ and by the product rule $(ab)'=a'b+ab'$. For $n=\prod p^{e_p}$, this implies

$$n'=n\sum_{p^{e_p}\parallel n}\frac{e_p}{p}.$$

The target value uses $N=5\times 10^{15}$, so factoring every integer up to $N$ is completely infeasible. The solution therefore turns $\gcd(n,n')$ into a multiplicative quantity and sums it by grouping integers according to their repeated prime factors.

Mathematical Approach

Write

$$g(n)=\gcd(n,n').$$

The whole problem is to compute $\sum_{n\le N} g(n)$ and then exclude $n=1$.

Step 1: Find the Local Prime-Power Factor

Fix a prime $p$, and write $n=p^e m$ with $p\nmid m$. Using the product rule,

$$n'=(p^e m)'=e\,p^{e-1}m+p^e m'=p^{e-1}(e\,m+p\,m').$$

Because $m$ is not divisible by $p$, the $p$-adic valuation of $n'$ behaves in a simple way:

$$v_p(n')=\begin{cases} e-1,& p\nmid e,\\ \ge e,& p\mid e. \end{cases}$$

Therefore the $p$-part of $\gcd(n,n')$ is

$$v_p(g(n))=\min\!\bigl(e,v_p(n')\bigr)=\begin{cases} e-1,& p\nmid e,\\ e,& p\mid e. \end{cases}$$

Equivalently, for one prime power,

$$g(p^e)=\begin{cases} p^{e-1},& p\nmid e,\\ p^e,& p\mid e. \end{cases}$$

Since each prime valuation is determined independently, $g$ is multiplicative:

$$g(n)=\prod_{p^e\parallel n} g(p^e).$$

Step 2: Separate the Powerful Part from the Squarefree Tail

If a prime appears only once, then $g(p)=1$. So primes of exponent $1$ do not change $g(n)$. This suggests writing every integer uniquely as

$$n=u\,s,$$

where $u$ is the product of all prime powers $p^e$ with $e\ge 2$, and $s$ is the product of the remaining primes of exponent $1$.

Then:

$$u \text{ is powerful},\qquad s \text{ is squarefree},\qquad \gcd(s,u)=1,$$

and, crucially,

$$g(n)=g(u).$$

If we define the radical of $u$ by

$$\operatorname{rad}(u)=\prod_{p\mid u} p,$$

then for each fixed powerful part $u$, the admissible tails are exactly the squarefree integers $s\le N/u$ that are coprime to $\operatorname{rad}(u)$. Hence

$$S(N)+1=\sum_{\substack{u\le N\\ u\ \mathrm{powerful}}} g(u)\,Q\!\left(\left\lfloor\frac{N}{u}\right\rfloor,\operatorname{rad}(u)\right),$$

where

$$Q(x,r)=\#\{\,s\le x:\ s\text{ is squarefree and }\gcd(s,r)=1\,\}.$$

The extra $1$ comes from $n=1$, which is handled separately at the end.

Step 3: Count the Squarefree Coprime Tail

The indicator of squarefreeness is

$$1_{\mathrm{sqf}}(s)=\sum_{d^2\mid s}\mu(d),$$

so Möbius inversion gives

$$Q(x,r)=\sum_{\substack{d\le \sqrt{x}\\ \gcd(d,r)=1}}\mu(d)\,\Phi_r\!\left(\left\lfloor\frac{x}{d^2}\right\rfloor\right),$$

where $\Phi_r(y)$ counts integers up to $y$ that are coprime to $r$:

$$\Phi_r(y)=\#\{\,m\le y:\gcd(m,r)=1\,\}.$$

Because $r=\operatorname{rad}(u)$ is squarefree, inclusion-exclusion over the primes dividing $r$ gives

$$\Phi_r(y)=\sum_{t\mid r}\mu(t)\left\lfloor\frac{y}{t}\right\rfloor.$$

This is the exact counting formula used by the full implementation: one outer Möbius sum for squarefree numbers, and one inner divisor sum for the coprimality condition.

Step 4: Enumerate All Powerful Parts Efficiently

Every powerful number can be built by choosing distinct primes and assigning each chosen prime an exponent at least $2$. A depth-first search over increasing primes therefore visits each powerful part $u\le N$ exactly once.

At each node, the algorithm knows three pieces of information:

$u$ itself, the already accumulated value $g(u)$, and the prime set of $\operatorname{rad}(u)$. With those data, it immediately adds

$$g(u)\,Q\!\left(\left\lfloor\frac{N}{u}\right\rfloor,\operatorname{rad}(u)\right)$$

to the answer, then extends $u$ by appending a new prime with exponent $2,3,4,\dots$ as long as the product stays within $N$.

The shorter implementations encode the same idea in an equivalent telescoping recursion. If $G_e(p)=g(p^e)$, then each increase

$$\Delta_e(p)=G_e(p)-G_{e-1}(p)$$

is multiplied by the number of integers whose powerful part first acquires the exponent $e$ at that prime, and the recursion continues with larger primes only. This compact recursion and the explicit $Q(x,r)$ formulation compute the same sum.

Worked Example: $N=20$

The powerful parts not exceeding $20$ are

$$1,\ 4,\ 8,\ 9,\ 16.$$

Now count the squarefree tails for each one:

$$\begin{aligned} u=1&:&&g(u)=1,\quad Q(20,1)=13,&&\text{contribution }13,\\ u=4&:&&g(u)=4,\quad Q(5,2)=3,&&\text{contribution }12,\\ u=8&:&&g(u)=4,\quad Q(2,2)=1,&&\text{contribution }4,\\ u=9&:&&g(u)=3,\quad Q(2,3)=2,&&\text{contribution }6,\\ u=16&:&&g(u)=16,\quad Q(1,2)=1,&&\text{contribution }16. \end{aligned}$$

Therefore

$$\sum_{n=1}^{20} g(n)=13+12+4+6+16=51,$$

so after removing the $n=1$ term we get

$$\sum_{n=2}^{20}\gcd(n,n')=50.$$

This small checkpoint matches the multiplicative formula perfectly.

How the Code Works

The C++, Python, and Java implementations all begin from the same prime-power formula for $g(p^e)$. The C++ implementation first validates that formula on small ranges, then sieves primes and Möbius values up to $\sqrt{N}$. It enumerates powerful parts by depth-first search and, for each one, counts squarefree coprime tails with the Möbius-plus-inclusion-exclusion formula above. The same file also contains a more compact equivalent recursion, and that is the version mirrored by the Python and Java implementations: they recurse over increasing primes, let the exponent at one prime grow from $2$ upward, and whenever the local factor of $g$ increases they add that increment times the number of admissible remaining multiples. All three languages therefore compute the same decomposition, just with different levels of explicit counting infrastructure.

Complexity Analysis

Let $M=\lfloor\sqrt{N}\rfloor$. The sieve stage costs $O(M\log\log M)$ time and $O(M)$ memory. The search space of powerful numbers is sparse: there are far fewer powerful integers up to $N$ than ordinary integers, so the recursion visits only a tiny fraction of the interval $[1,N]$. For each visited powerful part, the squarefree-coprime count sums over squarefree $d\le \sqrt{N/u}$, while the inner inclusion-exclusion runs over divisors of $\operatorname{rad}(u)$, whose prime support is small. In practice this is many orders of magnitude faster than iterating over every $n\le N$, which is why $N=5\times 10^{15}$ becomes feasible.

Footnotes and References

Problem page: https://projecteuler.net/problem=484
Arithmetic derivative: Wikipedia - Arithmetic derivative
Möbius function: Wikipedia - Möbius function
Squarefree integer: Wikipedia - Squarefree integer
Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle

Problem 484 source code

C++

#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <numeric>
#include <string>
#include <thread>
#include <vector>

using namespace std;

namespace {

using u128 = unsigned __int128;

struct Context {
    uint64_t N = 0;
    const vector<uint32_t>* primes = nullptr;
    const vector<int8_t>* mu = nullptr;
    const vector<uint32_t>* sqf_d = nullptr;
};

static uint64_t isqrt_u64(uint64_t x) {
    long double approx = sqrtl(static_cast<long double>(x));
    uint64_t r = static_cast<uint64_t>(approx);
    while ((r + 1) * (r + 1) <= x) ++r;
    while (r * r > x) --r;
    return r;
}

static vector<int8_t> mobius_sieve(int max_n, vector<uint32_t>& primes_out) {
    vector<int8_t> mu(max_n + 1, 0);
    vector<uint8_t> is_comp(max_n + 1, 0);
    mu[1] = 1;
    for (int i = 2; i <= max_n; ++i) {
        if (!is_comp[i]) {
            primes_out.push_back(static_cast<uint32_t>(i));
            mu[i] = -1;
        }
        for (uint32_t p : primes_out) {
            uint64_t v = static_cast<uint64_t>(i) * p;
            if (v > static_cast<uint64_t>(max_n)) break;
            is_comp[static_cast<size_t>(v)] = 1;
            if (i % static_cast<int>(p) == 0) {
                mu[static_cast<size_t>(v)] = 0;
                break;
            }
            mu[static_cast<size_t>(v)] = static_cast<int8_t>(-mu[i]);
        }
    }
    return mu;
}

// Counts squarefree s <= x with gcd(s, r)=1 where r is squarefree and encoded by its divisors/signs.
static uint64_t count_squarefree_coprime(uint64_t x,
                                         const vector<uint32_t>& primes_in_r,
                                         const vector<uint64_t>& divs,
                                         const vector<int8_t>& signs,
                                         const vector<int8_t>& mu,
                                         const vector<uint32_t>& sqf_d) {
    if (x == 0) return 0;
    uint64_t limit = isqrt_u64(x);
    auto end_it = upper_bound(sqf_d.begin(), sqf_d.end(), static_cast<uint32_t>(limit));
    const size_t prime_cnt = primes_in_r.size();

    int64_t res = 0;
    uint64_t last_y = numeric_limits<uint64_t>::max();
    int64_t last_phi = 0;

    for (auto it = sqf_d.begin(); it != end_it; ++it) {
        uint64_t d = *it;
        bool ok = true;
        if (prime_cnt == 1) {
            ok = (d % primes_in_r[0] != 0);
        } else if (prime_cnt == 2) {
            ok = (d % primes_in_r[0] != 0) && (d % primes_in_r[1] != 0);
        } else if (prime_cnt == 3) {
            ok = (d % primes_in_r[0] != 0) && (d % primes_in_r[1] != 0) &&
                 (d % primes_in_r[2] != 0);
        } else if (prime_cnt == 4) {
            ok = (d % primes_in_r[0] != 0) && (d % primes_in_r[1] != 0) &&
                 (d % primes_in_r[2] != 0) && (d % primes_in_r[3] != 0);
        } else {
            for (uint32_t p : primes_in_r) {
                if (d % p == 0) {
                    ok = false;
                    break;
                }
            }
        }
        if (!ok) continue;

        const int8_t mu_d = mu[static_cast<size_t>(d)];
        uint64_t y = x / (d * d);
        if (y != last_y) {
            int64_t phi = 0;
            for (size_t i = 0; i < divs.size(); ++i) {
                phi += signs[i] * static_cast<int64_t>(y / divs[i]);
            }
            last_y = y;
            last_phi = phi;
        }
        res += mu_d * last_phi;
    }
    return static_cast<uint64_t>(res);
}

// For prime power p^e (e>=2): g(p^e)=p^(e-1) unless p | e, then g=p^e.
static uint64_t g_prime_power(uint64_t p, int e, uint64_t p_pow) {
    if (static_cast<uint64_t>(e) % p == 0) return p_pow;
    return p_pow / p;
}

static void dfs(int start_idx,
                uint64_t curr_u,
                uint64_t curr_g,
                vector<uint32_t>& primes_in_r,
                vector<uint64_t>& divs,
                vector<int8_t>& signs,
                const Context& ctx,
                u128& acc) {
    // Each node corresponds to one powerful part; squarefree coprime factors are counted separately.
    uint64_t x = ctx.N / curr_u;
    uint64_t q = count_squarefree_coprime(x, primes_in_r, divs, signs, *ctx.mu, *ctx.sqf_d);
    acc += static_cast<u128>(curr_g) * q;

    const auto& primes = *ctx.primes;
    for (int i = start_idx; i < static_cast<int>(primes.size()); ++i) {
        uint64_t p = primes[i];
        uint64_t p2 = p * p;
        if (curr_u > ctx.N / p2) break;

        size_t primes_size = primes_in_r.size();
        size_t divs_size = divs.size();
        primes_in_r.push_back(static_cast<uint32_t>(p));
        for (size_t j = 0; j < divs_size; ++j) {
            divs.push_back(divs[j] * p);
            signs.push_back(static_cast<int8_t>(-signs[j]));
        }

        uint64_t p_pow = p2;
        int e = 2;
        while (curr_u <= ctx.N / p_pow) {
            uint64_t new_u = curr_u * p_pow;
            uint64_t g_p = g_prime_power(p, e, p_pow);
            uint64_t new_g = curr_g * g_p;
            dfs(i + 1, new_u, new_g, primes_in_r, divs, signs, ctx, acc);
            if (p_pow > ctx.N / p) break;
            p_pow *= p;
            ++e;
        }

        primes_in_r.resize(primes_size);
        divs.resize(divs_size);
        signs.resize(divs_size);
    }
}

static void worker(const Context& ctx, int prime_limit, atomic<int>& next_idx, u128& acc) {
    const auto& primes = *ctx.primes;
    while (true) {
        int idx = next_idx.fetch_add(1);
        if (idx >= prime_limit) break;
        uint64_t p = primes[idx];
        uint64_t p2 = p * p;
        if (p2 > ctx.N) continue;

        vector<uint32_t> primes_in_r;
        primes_in_r.reserve(8);
        primes_in_r.push_back(static_cast<uint32_t>(p));

        vector<uint64_t> divs;
        vector<int8_t> signs;
        divs.reserve(128);
        signs.reserve(128);
        divs.push_back(1);
        divs.push_back(p);
        signs.push_back(1);
        signs.push_back(-1);

        uint64_t p_pow = p2;
        int e = 2;
        while (p_pow <= ctx.N) {
            uint64_t g_p = g_prime_power(p, e, p_pow);
            dfs(idx + 1, p_pow, g_p, primes_in_r, divs, signs, ctx, acc);
            if (p_pow > ctx.N / p) break;
            p_pow *= p;
            ++e;
        }
    }
}

static u128 compute_sum(uint64_t N, bool use_threads) {
    uint64_t sqrt_n = isqrt_u64(N);
    vector<uint32_t> primes;
    vector<int8_t> mu = mobius_sieve(static_cast<int>(sqrt_n), primes);
    vector<uint32_t> sqf_d;
    sqf_d.reserve(static_cast<size_t>(sqrt_n * 7 / 10));
    for (uint32_t d = 1; d <= sqrt_n; ++d) {
        if (mu[static_cast<size_t>(d)] != 0) sqf_d.push_back(d);
    }

    Context ctx;
    ctx.N = N;
    ctx.primes = &primes;
    ctx.mu = &mu;
    ctx.sqf_d = &sqf_d;

    u128 total = 0;
    vector<uint32_t> primes_in_r;
    vector<uint64_t> divs = {1};
    vector<int8_t> signs = {1};
    uint64_t q_root = count_squarefree_coprime(N, primes_in_r, divs, signs, mu, sqf_d);
    total += q_root;

    int prime_limit = 0;
    while (prime_limit < static_cast<int>(primes.size())) {
        uint64_t p = primes[prime_limit];
        if (p * p > N) break;
        ++prime_limit;
    }

    int hw_threads = static_cast<int>(thread::hardware_concurrency());
    if (!use_threads || hw_threads <= 1 || prime_limit == 0) {
        atomic<int> next_idx(0);
        worker(ctx, prime_limit, next_idx, total);
    } else {
        int thread_count = hw_threads;
        atomic<int> next_idx(0);
        vector<thread> threads;
        vector<u128> partials(thread_count, 0);
        threads.reserve(thread_count);
        for (int t = 0; t < thread_count; ++t) {
            threads.emplace_back([&ctx, prime_limit, &next_idx, &partials, t]() {
                worker(ctx, prime_limit, next_idx, partials[t]);
            });
        }
        for (auto& th : threads) th.join();
        for (int t = 0; t < thread_count; ++t) {
            total += partials[t];
        }
    }
    return total;
}

static vector<uint32_t> generate_primes_upto(uint32_t limit) {
    vector<uint32_t> primes;
    vector<uint8_t> is_comp(static_cast<size_t>(limit + 1), 0);
    for (uint32_t i = 2; i <= limit; ++i) {
        if (!is_comp[static_cast<size_t>(i)]) {
            primes.push_back(i);
            if (static_cast<uint64_t>(i) * i <= limit) {
                for (uint64_t j = static_cast<uint64_t>(i) * i; j <= limit; j += i) {
                    is_comp[static_cast<size_t>(j)] = 1;
                }
            }
        }
    }
    return primes;
}

static u128 dfs_alt(int i0, uint64_t L0, const vector<uint32_t>& primes, const vector<uint64_t>& p2) {
    u128 res = 0;
    for (int i = i0; i < static_cast<int>(primes.size()); ++i) {
        const uint64_t q = p2[static_cast<size_t>(i)];
        uint64_t L = L0 / q;
        if (L == 0) break;

        const uint64_t p = primes[static_cast<size_t>(i)];
        uint64_t e = 1;
        uint64_t g = 1;
        while (L > 0) {
            const uint64_t gp = g;
            ++e;
            if (e != 1) {
                if (e == p) {
                    g *= q;
                    e = 0;
                } else {
                    g *= p;
                }
                const uint64_t c = g - gp;
                res += static_cast<u128>(c) * static_cast<u128>(L);
                if (L > q) {
                    res += static_cast<u128>(c) * dfs_alt(i + 1, L, primes, p2);
                }
            }
            L /= p;
        }
    }
    return res;
}

static u128 compute_sum_alt(uint64_t L) {
    const uint32_t limit = static_cast<uint32_t>(isqrt_u64(L));
    const vector<uint32_t> primes = generate_primes_upto(limit);
    vector<uint64_t> p2;
    p2.reserve(primes.size());
    for (uint32_t p : primes) p2.push_back(static_cast<uint64_t>(p) * p);
    return static_cast<u128>(L) + dfs_alt(0, L, primes, p2);
}

static vector<int> build_spf(int n) {
    vector<int> spf(n + 1, 0);
    for (int i = 2; i <= n; ++i) {
        if (spf[i] != 0) continue;
        spf[i] = i;
        if (static_cast<int64_t>(i) * i > n) continue;
        for (int j = i * i; j <= n; j += i) {
            if (spf[j] == 0) spf[j] = i;
        }
    }
    return spf;
}

static uint64_t arithmetic_derivative(uint64_t n, const vector<int>& spf) {
    if (n == 1) return 0;
    uint64_t x = n;
    uint64_t sum = 0;
    while (x > 1) {
        int p = spf[static_cast<size_t>(x)];
        if (p == 0) p = static_cast<int>(x);
        int e = 0;
        while (x % p == 0) {
            x /= p;
            ++e;
        }
        sum += static_cast<uint64_t>(e) * (n / p);
    }
    return sum;
}

static uint64_t gcd_formula_value(uint64_t n, const vector<int>& spf) {
    if (n == 1) return 1;
    uint64_t x = n;
    uint64_t g = 1;
    while (x > 1) {
        int p = spf[static_cast<size_t>(x)];
        if (p == 0) p = static_cast<int>(x);
        int e = 0;
        while (x % p == 0) {
            x /= p;
            ++e;
        }
        uint64_t p_pow = 1;
        for (int i = 1; i < e; ++i) p_pow *= static_cast<uint64_t>(p);
        uint64_t term = p_pow;
        if (static_cast<uint64_t>(e) % static_cast<uint64_t>(p) == 0) {
            term *= static_cast<uint64_t>(p);
        }
        g *= term;
    }
    return g;
}

static uint64_t brute_sum(int n) {
    vector<int> spf = build_spf(n);
    uint64_t sum = 0;
    for (int i = 1; i <= n; ++i) {
        uint64_t deriv = arithmetic_derivative(static_cast<uint64_t>(i), spf);
        sum += gcd<uint64_t>(static_cast<uint64_t>(i), deriv);
    }
    return sum;
}

static void validate_formula() {
    const int limit = 5000;
    vector<int> spf = build_spf(limit);
    for (int n = 1; n <= limit; ++n) {
        uint64_t deriv = arithmetic_derivative(static_cast<uint64_t>(n), spf);
        uint64_t g = gcd<uint64_t>(static_cast<uint64_t>(n), deriv);
        uint64_t g_formula = gcd_formula_value(static_cast<uint64_t>(n), spf);
        if (g != g_formula) {
            cerr << "Formula check failed at n=" << n << '\n';
            exit(1);
        }
    }
}

static void validate_small() {
    const int n = 10000;
    uint64_t brute = brute_sum(n);
    u128 fast = compute_sum(static_cast<uint64_t>(n), false);
    if (fast != static_cast<u128>(brute)) {
        cerr << "Small-N check failed: expected " << brute << '\n';
        exit(1);
    }
    u128 alt = compute_sum_alt(static_cast<uint64_t>(n));
    if (alt != static_cast<u128>(brute)) {
        cerr << "Small-N alt check failed: expected " << brute
             << " got " << static_cast<uint64_t>(alt) << '\n';
        exit(1);
    }
}

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

}  // namespace

int main() {
    validate_formula();
    validate_small();

    const uint64_t N = 5000000000000000ULL;
    u128 total = compute_sum_alt(N);
    if (total > 0) total -= 1;  // exclude n=1
    cout << to_string_u128(total) << '\n';
    return 0;
}

Python

import math

def solve():
    N = 5_000_000_000_000_000

    def isqrt(x): return math.isqrt(x)

    limit = isqrt(N)
    # Sieve primes up to sqrt(N)
    is_comp = bytearray(limit + 1)
    primes = []
    for i in range(2, limit + 1):
        if not is_comp[i]:
            primes.append(i)
            if i <= limit // i:
                for j in range(i*i, limit+1, i):
                    is_comp[j] = 1
    p2 = [p * p for p in primes]

    def dfs_alt(i0, L0):
        res = 0
        for i in range(i0, len(primes)):
            q = p2[i]
            L = L0 // q
            if L == 0: break
            p = primes[i]
            e = 1
            g = 1
            while L > 0:
                gp = g
                e += 1
                if e != 1:
                    if e == p:
                        g *= q
                        e = 0
                    else:
                        g *= p
                    c = g - gp
                    res += c * L
                    if L > q:
                        res += c * dfs_alt(i + 1, L)
                L //= p
        return res

    total = N + dfs_alt(0, N)
    total -= 1  # exclude n=1
    return str(total)

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

Java

import java.util.ArrayList;
import java.util.List;

public class Euler484 {

    private static List<Long> primes = new ArrayList<>();
    private static List<Long> p2 = new ArrayList<>();

    private static void generatePrimesUpto(int limit) {
        byte[] isComp = new byte[limit + 1];
        for (int i = 2; i <= limit; ++i) {
            if (isComp[i] == 0) {
                primes.add((long) i);
                p2.add((long) i * i);
                if ((long) i * i <= limit) {
                    for (long j = (long) i * i; j <= limit; j += i) {
                        isComp[(int) j] = 1;
                    }
                }
            }
        }
    }

    private static long dfsAlt(int i0, long L0) {
        long res = 0;
        for (int i = i0; i < primes.size(); ++i) {
            long q = p2.get(i);
            long L = L0 / q;
            if (L == 0)
                break;

            long p = primes.get(i);
            long e = 1;
            long g = 1;
            while (L > 0) {
                long gp = g;
                e++;
                if (e != 1) {
                    if (e == p) {
                        g *= q;
                        e = 0;
                    } else {
                        g *= p;
                    }
                    long c = g - gp;
                    res += c * L;
                    if (L > q) {
                        res += c * dfsAlt(i + 1, L);
                    }
                }
                L /= p;
            }
        }
        return res;
    }

    private static long computeSumAlt(long L) {
        int limit = (int) Math.sqrt(L);
        generatePrimesUpto(limit);
        return L + dfsAlt(0, L);
    }

    public static void main(String[] args) {
        long N = 5000000000000000L;
        long total = computeSumAlt(N);
        if (total > 0) {
            total -= 1;
        }
        System.out.println(total);
    }
}