Problem 850: Fractions of Powers

Project Euler Problem 850 Solution

EulerSolve provides an optimized solution for Project Euler Problem 850, Fractions of Powers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For odd $k$, define $$f_k(n)=\sum_{i=1}^{n}\left\{\frac{i^k}{n}\right\},$$ where $\{x\}$ denotes the fractional part of $x$. The problem then asks for $$S(N)=\sum_{\substack{k=1\\k\text{ odd}}}^{N}\sum_{n=1}^{N} f_k(n),$$ and specifically for $\lfloor S(33557799775533)\rfloor \bmod 977676779$. A direct double loop over all odd $k$ and all $n$ is hopeless at this scale, so the solution turns the fractional-part sum into a much more structured divisibility problem. Mathematical Approach The C++, Python, and Java implementations all follow the same chain of reductions: first remove the fractional parts, then express the remaining term through prime exponents, and finally reorganize the whole computation by divisor profiles. Step 1: Pair $i$ with $n-i$ Because $k$ is odd, $$ (n-i)^k\equiv -\,i^k \pmod n. $$ So for $1\le i\le n$, the two terms $\left\{\frac{i^k}{n}\right\}$ and $\left\{\frac{(n-i)^k}{n}\right\}$ add up to $1$, except when $n\mid i^k$, in which case both fractional parts are $0$. Define $$z_k(n)=\#\left\{1\le i\le n:\ n\mid i^k\right\}.$$ Then the pairing identity becomes $$2f_k(n)=n-z_k(n).$$ This is the first crucial simplification: the problem is no longer about fractional parts themselves, only about how many residues produce exact divisibility....

Detailed mathematical approach

Problem Summary

For odd $k$, define

$$f_k(n)=\sum_{i=1}^{n}\left\{\frac{i^k}{n}\right\},$$

where $\{x\}$ denotes the fractional part of $x$. The problem then asks for

$$S(N)=\sum_{\substack{k=1\\k\text{ odd}}}^{N}\sum_{n=1}^{N} f_k(n),$$

and specifically for $\lfloor S(33557799775533)\rfloor \bmod 977676779$. A direct double loop over all odd $k$ and all $n$ is hopeless at this scale, so the solution turns the fractional-part sum into a much more structured divisibility problem.

Mathematical Approach

The C++, Python, and Java implementations all follow the same chain of reductions: first remove the fractional parts, then express the remaining term through prime exponents, and finally reorganize the whole computation by divisor profiles.

Step 1: Pair $i$ with $n-i$

Because $k$ is odd,

$$ (n-i)^k\equiv -\,i^k \pmod n. $$

So for $1\le i\le n$, the two terms $\left\{\frac{i^k}{n}\right\}$ and $\left\{\frac{(n-i)^k}{n}\right\}$ add up to $1$, except when $n\mid i^k$, in which case both fractional parts are $0$. Define

$$z_k(n)=\#\left\{1\le i\le n:\ n\mid i^k\right\}.$$

Then the pairing identity becomes

$$2f_k(n)=n-z_k(n).$$

This is the first crucial simplification: the problem is no longer about fractional parts themselves, only about how many residues produce exact divisibility.

Step 2: Count the divisible residues prime by prime

Write

$$n=\prod_{p} p^{e_p}.$$

The condition $n\mid i^k$ is equivalent to

$$v_p(i)\ge \left\lceil\frac{e_p}{k}\right\rceil\qquad\text{for every prime }p\mid n,$$

where $v_p$ is the $p$-adic valuation. Therefore $i$ must be a multiple of

$$m_k(n)=\prod_{p\mid n} p^{\left\lceil e_p/k\right\rceil}.$$

Among the integers $1,2,\dots,n$, exactly $n/m_k(n)$ are such multiples. Hence

$$z_k(n)=\frac{n}{m_k(n)}$$

and therefore

$$2f_k(n)=n-\frac{n}{m_k(n)}.$$

Step 3: Split off the easy closed form

Define

$$T(N)=2S(N).$$

Using the previous identity,

$$T(N)=\sum_{\substack{k=1\\k\text{ odd}}}^{N}\sum_{n=1}^{N}\left(n-\frac{n}{m_k(n)}\right)=A(N)-H(N),$$

where

$$A(N)=\sum_{\substack{k=1\\k\text{ odd}}}^{N}\sum_{n=1}^{N} n=\frac{N+1}{2}\cdot\frac{N(N+1)}{2}$$

for the actual odd input $N$, and

$$H(N)=\sum_{\substack{k=1\\k\text{ odd}}}^{N}\sum_{n=1}^{N}\frac{n}{m_k(n)}.$$

So the entire problem has been reduced to evaluating $H(N)$ efficiently.

Step 4: Rewrite $\frac{n}{m_k(n)}$ as a totient-weighted divisor sum

The implementation does not sum $\frac{n}{m_k(n)}$ directly. Instead, for odd $k\ge 3$, it introduces

$$\rho_k(d)=\prod_{p^a\parallel d} p^{\left\lceil a/(k-1)\right\rceil}.$$

Then one has the identity

$$\frac{n}{m_k(n)}=\sum_{\substack{d\ge 1\\ d\,\rho_k(d)\mid n}} \varphi(d),$$

where $\varphi$ is Euler's totient function.

Why is this true? It is enough to check one prime power. If $n=p^E$, then the admissible exponents $a$ in $d=p^a$ are exactly those satisfying

$$a+\left\lceil\frac{a}{k-1}\right\rceil\le E.$$

The largest such $a$ is $E-\left\lceil E/k\right\rceil$, so

$$\sum_{\substack{a\ge 0\\ a+\lceil a/(k-1)\rceil\le E}} \varphi(p^a)=\sum_{a=0}^{E-\lceil E/k\rceil}\varphi(p^a)=p^{E-\lceil E/k\rceil}=\frac{p^E}{p^{\lceil E/k\rceil}}.$$

Since both sides are multiplicative in $n$, the formula follows for general $n$. For $k=1$, the situation is simpler: $m_1(n)=n$, hence $\frac{n}{m_1(n)}=1$.

Step 5: Compress all odd $k$ with the same ceiling pattern

Substituting the divisor identity gives

$$H(N)=N+\sum_{\substack{k=3\\k\text{ odd}}}^{N}\ \sum_{d\ge 1}\varphi(d)\left\lfloor\frac{N}{d\,\rho_k(d)}\right\rfloor.$$

Now fix one divisor profile

$$d=\prod_{j=1}^{r} p_j^{a_j}.$$

Only profiles with

$$d\,\operatorname{rad}(d)\le N,\qquad \operatorname{rad}(d)=\prod_{p\mid d}p,$$

can contribute, because $\rho_k(d)\ge \operatorname{rad}(d)$ for every odd $k\ge 3$. This immediately implies that every prime in $d$ is at most $\sqrt N$, which is why the implementations only sieve primes up to $\sqrt N$.

For a fixed $d$, the values $\left\lceil a_j/(k-1)\right\rceil$ do not change at every odd $k$; they change only when $k-1$ crosses one of finitely many exponent thresholds. Once $k$ is larger than the biggest exponent $a_j$, every ceiling becomes $1$, so

$$\rho_k(d)=\operatorname{rad}(d).$$

That means the whole tail of large odd $k$ can be summed in one block, while the small odd $k$ values are handled individually. This is the compression that makes the problem tractable.

Worked Example: $n=72$ and $k=3$

Take

$$72=2^3\cdot 3^2.$$

For $k=3$,

$$m_3(72)=2^{\lceil 3/3\rceil}3^{\lceil 2/3\rceil}=2\cdot 3=6,$$

$$z_3(72)=\frac{72}{6}=12,\qquad 2f_3(72)=72-12=60,\qquad f_3(72)=30.$$

The divisor-profile formula gives the same result. Here

$$\rho_3(d)=\prod_{p^a\parallel d} p^{\lceil a/2\rceil}.$$

For the prime $2$, the admissible exponents $a$ satisfy $a+\lceil a/2\rceil\le 3$, so $a\in\{0,1,2\}$. For the prime $3$, the condition is $a+\lceil a/2\rceil\le 2$, so $a\in\{0,1\}$. Therefore

$$\sum_{\substack{d\,\rho_3(d)\mid 72}}\varphi(d)=\left(1+\varphi(2)+\varphi(4)\right)\left(1+\varphi(3)\right)=(1+1+2)(1+2)=12,$$

which matches $\frac{72}{6}$ exactly.

How the Code Works

The implementations first generate all primes up to $\sqrt N$ with a linear sieve. That is enough because every contributing divisor profile $d$ must satisfy $d\,\operatorname{rad}(d)\le N$, so no prime factor larger than $\sqrt N$ can appear in a nontrivial profile.

They then enumerate divisor profiles $d=\prod p^a$ recursively in increasing prime order. During that recursion they maintain four multiplicative pieces of state: the value of $d$, the radical $\operatorname{rad}(d)$, the totient $\varphi(d)$, and the largest exponent currently present in the profile. The recursion stops immediately once $d\,\operatorname{rad}(d)>N$, which is the main pruning rule.

For each profile, the code evaluates

$$\sum_{\substack{k=3\\k\text{ odd}}}^{N}\left\lfloor\frac{N}{d\,\rho_k(d)}\right\rfloor$$

without iterating over every odd $k$. Large odd exponents share the same denominator $d\,\operatorname{rad}(d)$, so they are aggregated in one arithmetic block; only the finitely many smaller breakpoints are processed one by one. After multiplying by $\varphi(d)$ and summing all profiles, the implementation obtains $H(N)$, computes $T(N)=A(N)-H(N)$ modulo $2M$ with $M=977676779$, and finally converts that residue into

$$\left\lfloor\frac{T(N)}{2}\right\rfloor \bmod M=\lfloor S(N)\rfloor \bmod M.$$

The C++ and Java implementations also split the top-level prime branches across threads, while the Python implementation uses the same arithmetic in a single-threaded recursive traversal.

Complexity Analysis

The prime sieve up to $\sqrt N$ costs $O(\sqrt N)$ time and $O(\sqrt N)$ memory. The harder part is the recursive enumeration of divisor profiles. If

$$\mathcal{D}(N)=\left\{d\ge 1:\ d\,\operatorname{rad}(d)\le N\right\},$$

then the recursion visits exactly the admissible profiles in $\mathcal{D}(N)$, and for each profile it evaluates only the finitely many odd-$k$ breakpoint intervals where some ceiling $\left\lceil a/(k-1)\right\rceil$ changes. So the running time is roughly

$$O\!\left(\sqrt N+\sum_{d\in\mathcal{D}(N)}(1+B(d))\right),$$

where $B(d)$ is the number of distinct odd-$k$ breakpoints induced by the exponents of $d$. There is no simple closed elementary form for this quantity, but it is dramatically smaller than the original $O(N^2)$ scan over all $(k,n)$. Memory usage stays at $O(\sqrt N)$ for the prime table plus $O(\log N)$ recursion depth.

Footnotes and References

Problem page: https://projecteuler.net/problem=850
Fractional part: Wikipedia — Fractional part
$p$-adic valuation: Wikipedia — p-adic valuation
Euler's totient function: Wikipedia — Euler's totient function
Ceiling function: Wikipedia — Floor and ceiling functions

Problem 850 source code

C++

#include <algorithm>
#include <atomic>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <pthread.h>
#include <unistd.h>
#include <vector>

using u32 = uint32_t;
using u64 = uint64_t;
using i64 = int64_t;
using u128 = unsigned __int128;
using i128 = __int128_t;

namespace {

constexpr u64 MOD = 977676779ULL;

u64 add_mod(u64 a, u64 b, u64 mod) {
    a += b;
    if (a >= mod) {
        a -= mod;
    }
    return a;
}

u64 mul_mod(u64 a, u64 b, u64 mod) {
    return static_cast<u64>((static_cast<u128>(a) * b) % mod);
}

u64 norm_i128_mod(i128 x, u64 mod) {
    i128 r = x % static_cast<i128>(mod);
    if (r < 0) {
        r += static_cast<i128>(mod);
    }
    return static_cast<u64>(r);
}

u64 isqrt_u64(u64 x) {
    u64 lo = 0;
    u64 hi = (x < (1ULL << 32)) ? (1ULL << 32) : x;
    while (lo + 1 < hi) {
        u64 mid = lo + (hi - lo) / 2;
        if (mid <= x / mid) {
            lo = mid;
        } else {
            hi = mid;
        }
    }
    return lo;
}

std::vector<u32> primes_up_to(int n) {
    std::vector<int> lp(n + 1, 0);
    std::vector<u32> primes;
    primes.reserve(static_cast<size_t>(n / 10));
    for (int i = 2; i <= n; ++i) {
        if (lp[i] == 0) {
            lp[i] = i;
            primes.push_back(static_cast<u32>(i));
        }
        for (u32 p : primes) {
            i64 v = static_cast<i64>(i) * static_cast<i64>(p);
            if (v > n || static_cast<int>(p) > lp[i]) {
                break;
            }
            lp[static_cast<size_t>(v)] = static_cast<int>(p);
        }
    }
    return primes;
}

u64 pow_int(u64 base, int exp) {
    u64 res = 1;
    for (int i = 0; i < exp; ++i) {
        res *= base;
    }
    return res;
}

u64 brute_T_small(int N) {
    std::vector<int> spf(N + 1, 0);
    for (int i = 2; i <= N; ++i) {
        if (spf[i] == 0) {
            spf[i] = i;
            if (static_cast<i64>(i) * i <= N) {
                for (int j = i * i; j <= N; j += i) {
                    if (spf[j] == 0) {
                        spf[j] = i;
                    }
                }
            }
        }
    }

    std::vector<std::vector<std::pair<int, int>>> factors(N + 1);
    for (int n = 2; n <= N; ++n) {
        int x = n;
        while (x > 1) {
            int p = spf[x];
            int e = 0;
            while (x % p == 0) {
                x /= p;
                ++e;
            }
            factors[n].push_back({p, e});
        }
    }

    u64 T = 0;
    for (int k = 1; k <= N; k += 2) {
        for (int n = 1; n <= N; ++n) {
            u64 m = 1;
            for (const auto& pe : factors[n]) {
                int c = (pe.second + k - 1) / k;
                m *= pow_int(static_cast<u64>(pe.first), c);
            }
            T += static_cast<u64>(n - n / static_cast<int>(m));
        }
    }
    return T;
}

class Solver850 {
public:
    Solver850(u64 n, u64 mod, bool use_threads)
        : N_(n),
          mod2_(2 * mod),
          odd_half_((n + 1) / 2),
          use_threads_(use_threads) {
        const int lim = static_cast<int>(isqrt_u64(N_)) + 1;
        primes_ = primes_up_to(lim);
        prime_limit_ = 0;
        while (prime_limit_ < static_cast<int>(primes_.size()) &&
               static_cast<u64>(primes_[prime_limit_]) * static_cast<u64>(primes_[prime_limit_]) <= N_) {
            ++prime_limit_;
        }
    }

    u64 compute_T_mod() const {
        std::vector<u32> ps;
        std::vector<int> es;
        u64 h_mod = contribute(1ULL, 1ULL, ps, es, -1, 1ULL);
        h_mod = add_mod(h_mod, compute_branches_mod(), mod2_);

        const u64 sum_n_mod = static_cast<u64>((static_cast<u128>(N_) * (N_ + 1) / 2) % mod2_);
        const u64 a_mod = mul_mod(odd_half_ % mod2_, sum_n_mod, mod2_);
        if (a_mod >= h_mod) {
            return a_mod - h_mod;
        }
        return a_mod + mod2_ - h_mod;
    }

private:
    u64 contribute(u64 prd,
                   u64 rad,
                   const std::vector<u32>& ps,
                   const std::vector<int>& es,
                   int max_e,
                   u64 phi_mod) const {
        const u64 Np = N_ / prd;
        const int K = max_e + 2;
        const i64 shift = static_cast<i64>(K / 2);
        const i128 base = static_cast<i128>(Np / rad) * static_cast<i128>(static_cast<i64>(odd_half_) - shift);
        u64 ret = norm_i128_mod(base, mod2_);

        for (int k = K - (K & 1) - 2; k > 0; k -= 2) {
            u64 newr = 1;
            bool over = false;
            for (size_t i = 0; i < ps.size(); ++i) {
                const u64 p = ps[i];
                const int need = es[i] / k + 1;
                for (int t = 0; t < need; ++t) {
                    if (newr > Np / p) {
                        over = true;
                        break;
                    }
                    newr *= p;
                }
                if (over) {
                    break;
                }
            }
            if (over) {
                break;
            }
            ret += Np / newr;
            ret %= mod2_;
        }

        return mul_mod(ret, phi_mod, mod2_);
    }

    u64 rec(int ind,
            u64 prd,
            u64 rad,
            std::vector<u32>& ps,
            std::vector<int>& es,
            int max_e,
            u64 phi_mod) const {
        u64 ret = contribute(prd, rad, ps, es, max_e, phi_mod);

        for (int i = ind; i < prime_limit_; ++i) {
            const u64 p = primes_[i];
            if (prd > N_ / p || rad > N_ / p) {
                break;
            }
            const u64 nxt_prd = prd * p;
            const u64 nxt_rad = rad * p;
            const u64 prd_lim = N_ / nxt_rad;
            if (nxt_prd > prd_lim) {
                break;
            }

            ps.push_back(static_cast<u32>(p));
            es.push_back(-1);

            u64 cur_prd = nxt_prd;
            u64 phi_e_mod = phi_mod;
            while (cur_prd <= prd_lim) {
                const int e = es.back() + 1;
                es.back() = e;
                if (e == 0) {
                    phi_e_mod = mul_mod(phi_mod, (p - 1) % mod2_, mod2_);
                } else {
                    phi_e_mod = mul_mod(phi_e_mod, p % mod2_, mod2_);
                }

                const int child_max = std::max(max_e, e);
                ret = add_mod(ret, rec(i + 1, cur_prd, nxt_rad, ps, es, child_max, phi_e_mod), mod2_);

                if (cur_prd > prd_lim / p) {
                    break;
                }
                cur_prd *= p;
            }

            ps.pop_back();
            es.pop_back();
        }

        return ret;
    }

    u64 branch_from_prime(int idx) const {
        const u64 p = primes_[idx];
        if (p > N_ / p) {
            return 0;
        }

        const u64 prd_lim = N_ / p;
        if (p > prd_lim) {
            return 0;
        }

        std::vector<u32> ps(1, static_cast<u32>(p));
        std::vector<int> es(1, -1);
        u64 ret = 0;
        u64 cur_prd = p;
        u64 phi_e_mod = 1;
        while (cur_prd <= prd_lim) {
            const int e = es[0] + 1;
            es[0] = e;
            if (e == 0) {
                phi_e_mod = (p - 1) % mod2_;
            } else {
                phi_e_mod = mul_mod(phi_e_mod, p % mod2_, mod2_);
            }
            ret = add_mod(ret, rec(idx + 1, cur_prd, p, ps, es, e, phi_e_mod), mod2_);
            if (cur_prd > prd_lim / p) {
                break;
            }
            cur_prd *= p;
        }
        return ret;
    }

    struct WorkerArg {
        const Solver850* self = nullptr;
        std::atomic<int>* next_idx = nullptr;
        u64 partial = 0;
    };

    static void* worker_entry(void* raw) {
        auto* arg = static_cast<WorkerArg*>(raw);
        u64 local = 0;
        while (true) {
            const int idx = arg->next_idx->fetch_add(1, std::memory_order_relaxed);
            if (idx >= arg->self->prime_limit_) {
                break;
            }
            local = add_mod(local, arg->self->branch_from_prime(idx), arg->self->mod2_);
        }
        arg->partial = local;
        return nullptr;
    }

    u64 compute_branches_mod() const {
        if (prime_limit_ == 0) {
            return 0;
        }

        long cpu = ::sysconf(_SC_NPROCESSORS_ONLN);
        int thread_count = (cpu > 1) ? static_cast<int>(cpu) : 1;
        if (!use_threads_ || thread_count <= 1) {
            u64 total = 0;
            for (int i = 0; i < prime_limit_; ++i) {
                total = add_mod(total, branch_from_prime(i), mod2_);
            }
            return total;
        }

        if (thread_count > 8) {
            thread_count = 8;
        }
        if (thread_count > prime_limit_) {
            thread_count = prime_limit_;
        }

        std::atomic<int> next_idx(0);
        std::vector<pthread_t> tids(static_cast<size_t>(thread_count));
        std::vector<WorkerArg> args(static_cast<size_t>(thread_count));

        for (int t = 0; t < thread_count; ++t) {
            args[static_cast<size_t>(t)].self = this;
            args[static_cast<size_t>(t)].next_idx = &next_idx;
            const int rc = pthread_create(&tids[static_cast<size_t>(t)], nullptr, worker_entry, &args[static_cast<size_t>(t)]);
            assert(rc == 0);
        }

        u64 total = 0;
        for (int t = 0; t < thread_count; ++t) {
            const int rc = pthread_join(tids[static_cast<size_t>(t)], nullptr);
            assert(rc == 0);
            total = add_mod(total, args[static_cast<size_t>(t)].partial, mod2_);
        }
        return total;
    }

    u64 N_;
    u64 mod2_;
    u64 odd_half_;
    bool use_threads_;
    std::vector<u32> primes_;
    int prime_limit_;
};

void validate() {
    assert(brute_T_small(10) == 201ULL);
    assert(brute_T_small(30) == Solver850(30ULL, MOD, false).compute_T_mod());
    assert(Solver850(10ULL, MOD, false).compute_T_mod() == 201ULL);
    assert(Solver850(1000ULL, MOD, false).compute_T_mod() == 247375608ULL);
}

}  // namespace

int main() {
    validate();

    constexpr u64 N = 33557799775533ULL;
    const u64 t_mod = Solver850(N, MOD, true).compute_T_mod();
    const u64 answer = (t_mod - (t_mod & 1ULL)) / 2ULL;
    std::cout << (answer % MOD) << '\n';
    return 0;
}

Python

import math

def solve():
    MOD = 977676779
    N = 33557799775533
    mod2 = 2 * MOD; oh = (N+1)//2

    lim = int(N**0.5) + 1
    # Linear sieve for primes up to sqrt(N)
    lp = [0]*(lim+1); primes = []
    for i in range(2, lim+1):
        if lp[i] == 0: lp[i] = i; primes.append(i)
        for p in primes:
            if i*p > lim or p > lp[i]: break
            lp[i*p] = p
    pl = len(primes); prime_limit = 0
    while prime_limit < pl and primes[prime_limit]**2 <= N: prime_limit += 1

    def add_mod(a, b): return a+b-mod2 if a+b >= mod2 else a+b
    def mul_mod(a, b): return a*b % mod2

    def contribute(prd, rad, ps, es, max_e, phi_mod):
        Np = N // prd; K = max_e + 2; shift = K // 2
        base = (Np // rad) * (oh - shift)
        ret = base % mod2
        if ret < 0: ret += mod2
        for k in range(K - (K & 1) - 2, 0, -2):
            newr = 1; over = False
            for i in range(len(ps)):
                need = es[i] // k + 1
                for _ in range(need):
                    if newr > Np // ps[i]: over = True; break
                    newr *= ps[i]
                if over: break
            if over: break
            ret = (ret + Np // newr) % mod2
        return mul_mod(ret, phi_mod)

    def rec(ind, prd, rad, ps, es, max_e, phi_mod):
        ret = contribute(prd, rad, ps, es, max_e, phi_mod)
        for i in range(ind, prime_limit):
            p = primes[i]
            if prd > N // p or rad > N // p: break
            nxt_prd = prd * p; nxt_rad = rad * p
            prd_lim = N // nxt_rad
            if nxt_prd > prd_lim: break
            ps.append(p); es.append(-1)
            cur_prd = nxt_prd; phi_e = phi_mod
            while cur_prd <= prd_lim:
                e = es[-1] + 1; es[-1] = e
                if e == 0: phi_e = mul_mod(phi_mod, (p-1) % mod2)
                else: phi_e = mul_mod(phi_e, p % mod2)
                cm = max(max_e, e)
                ret = add_mod(ret, rec(i+1, cur_prd, nxt_rad, ps, es, cm, phi_e))
                if cur_prd > prd_lim // p: break
                cur_prd *= p
            ps.pop(); es.pop()
        return ret

    h_mod = contribute(1, 1, [], [], -1, 1)
    for i in range(prime_limit):
        p = primes[i]
        if p > N // p: continue
        prd_lim = N // p
        if p > prd_lim: continue
        ps = [p]; es = [-1]; cur_prd = p; phi_e = 1
        while cur_prd <= prd_lim:
            e = es[0] + 1; es[0] = e
            if e == 0: phi_e = (p-1) % mod2
            else: phi_e = mul_mod(phi_e, p % mod2)
            h_mod = add_mod(h_mod, rec(i+1, cur_prd, p, ps, es, e, phi_e))
            if cur_prd > prd_lim // p: break
            cur_prd *= p

    sum_n = N * (N+1) // 2 % mod2
    a_mod = (oh % mod2) * sum_n % mod2
    t_mod = (a_mod - h_mod + mod2) % mod2
    answer = (t_mod - (t_mod & 1)) // 2
    return str(answer % MOD)

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

Java

import java.util.ArrayList;
import java.util.concurrent.atomic.AtomicInteger;

public class Euler850 {
    static final long MOD = 977676779L;

    static long addMod(long a, long b, long mod) {
        a += b;
        if (a >= mod)
            a -= mod;
        return a;
    }

    static long mulMod(long a, long b, long mod) {
        java.math.BigInteger biA = java.math.BigInteger.valueOf(a);
        java.math.BigInteger biB = java.math.BigInteger.valueOf(b);
        java.math.BigInteger biMod = java.math.BigInteger.valueOf(mod);
        return biA.multiply(biB).remainder(biMod).longValue();
    }

    static long normI128Mod(java.math.BigInteger x, long mod) {
        java.math.BigInteger biMod = java.math.BigInteger.valueOf(mod);
        java.math.BigInteger r = x.remainder(biMod);
        if (r.compareTo(java.math.BigInteger.ZERO) < 0) {
            r = r.add(biMod);
        }
        return r.longValue();
    }

    static long isqrt(long x) {
        long lo = 0;
        long hi = (x < (1L << 32)) ? (1L << 32) : x;
        while (lo + 1 < hi) {
            long mid = lo + (hi - lo) / 2;
            if (mid <= x / mid) {
                lo = mid;
            } else {
                hi = mid;
            }
        }
        return lo;
    }

    static ArrayList<Long> primesUpTo(int n) {
        int[] lp = new int[n + 1];
        ArrayList<Long> primes = new ArrayList<>();
        for (int i = 2; i <= n; ++i) {
            if (lp[i] == 0) {
                lp[i] = i;
                primes.add((long) i);
            }
            for (long p : primes) {
                if (p * i > n || p > lp[i])
                    break;
                lp[(int) (p * i)] = (int) p;
            }
        }
        return primes;
    }

    static class Solver850 {
        long N;
        long mod2;
        long oddHalf;
        ArrayList<Long> primes;
        int primeLimit;

        Solver850(long n, long mod) {
            this.N = n;
            this.mod2 = 2 * mod;
            this.oddHalf = (n + 1) / 2;

            int lim = (int) isqrt(N) + 1;
            this.primes = primesUpTo(lim);
            this.primeLimit = 0;
            while (primeLimit < primes.size() && primes.get(primeLimit) * primes.get(primeLimit) <= N) {
                primeLimit++;
            }
        }

        long computeTMod() {
            ArrayList<Long> ps = new ArrayList<>();
            ArrayList<Integer> es = new ArrayList<>();
            long hMod = contribute(1L, 1L, ps, es, -1, 1L);
            hMod = addMod(hMod, computeBranchesMod(), mod2);

            java.math.BigInteger bn = java.math.BigInteger.valueOf(N);
            java.math.BigInteger sumNModBase = bn.multiply(bn.add(java.math.BigInteger.ONE))
                    .divide(java.math.BigInteger.valueOf(2));
            long sumNMod = sumNModBase.remainder(java.math.BigInteger.valueOf(mod2)).longValue();
            long aMod = mulMod(oddHalf % mod2, sumNMod, mod2);

            if (aMod >= hMod)
                return aMod - hMod;
            return aMod + mod2 - hMod;
        }

        long contribute(long prd, long rad, ArrayList<Long> ps, ArrayList<Integer> es, int maxE, long phiMod) {
            long Np = N / prd;
            int K = maxE + 2;
            long shift = K / 2;

            java.math.BigInteger base1 = java.math.BigInteger.valueOf(Np / rad);
            java.math.BigInteger base2 = java.math.BigInteger.valueOf(oddHalf - shift);
            java.math.BigInteger base = base1.multiply(base2);
            long ret = normI128Mod(base, mod2);

            for (int k = K - (K & 1) - 2; k > 0; k -= 2) {
                long newr = 1;
                boolean over = false;
                for (int i = 0; i < ps.size(); ++i) {
                    long p = ps.get(i);
                    int need = es.get(i) / k + 1;
                    for (int t = 0; t < need; ++t) {
                        if (newr > Np / p) {
                            over = true;
                            break;
                        }
                        newr *= p;
                    }
                    if (over)
                        break;
                }
                if (over)
                    break;

                ret += Np / newr;
                ret %= mod2;
            }

            return mulMod(ret, phiMod, mod2);
        }

        long rec(int ind, long prd, long rad, ArrayList<Long> ps, ArrayList<Integer> es, int maxE, long phiMod) {
            long ret = contribute(prd, rad, ps, es, maxE, phiMod);

            for (int i = ind; i < primeLimit; ++i) {
                long p = primes.get(i);
                if (prd > N / p || rad > N / p)
                    break;
                long nxtPrd = prd * p;
                long nxtRad = rad * p;
                long prdLim = N / nxtRad;
                if (nxtPrd > prdLim)
                    break;

                ps.add(p);
                es.add(-1);

                long curPrd = nxtPrd;
                long phiEMod = phiMod;
                while (curPrd <= prdLim) {
                    int e = es.get(es.size() - 1) + 1;
                    es.set(es.size() - 1, e);
                    if (e == 0) {
                        phiEMod = mulMod(phiMod, (p - 1) % mod2, mod2);
                    } else {
                        phiEMod = mulMod(phiEMod, p % mod2, mod2);
                    }

                    int childMax = Math.max(maxE, e);
                    ret = addMod(ret, rec(i + 1, curPrd, nxtRad, ps, es, childMax, phiEMod), mod2);

                    if (curPrd > prdLim / p)
                        break;
                    curPrd *= p;
                }

                ps.remove(ps.size() - 1);
                es.remove(es.size() - 1);
            }

            return ret;
        }

        long branchFromPrime(int idx) {
            long p = primes.get(idx);
            if (p > N / p)
                return 0;

            long prdLim = N / p;
            if (p > prdLim)
                return 0;

            ArrayList<Long> ps = new ArrayList<>();
            ps.add(p);
            ArrayList<Integer> es = new ArrayList<>();
            es.add(-1);

            long ret = 0;
            long curPrd = p;
            long phiEMod = 1;

            while (curPrd <= prdLim) {
                int e = es.get(0) + 1;
                es.set(0, e);
                if (e == 0) {
                    phiEMod = (p - 1) % mod2;
                } else {
                    phiEMod = mulMod(phiEMod, p % mod2, mod2);
                }

                ret = addMod(ret, rec(idx + 1, curPrd, p, ps, es, e, phiEMod), mod2);

                if (curPrd > prdLim / p)
                    break;
                curPrd *= p;
            }

            return ret;
        }

        long computeBranchesMod() {
            int cores = Runtime.getRuntime().availableProcessors();
            Thread[] threads = new Thread[cores];
            long[] partials = new long[cores];
            AtomicInteger nextIdx = new AtomicInteger(0);

            for (int t = 0; t < cores; ++t) {
                final int tid = t;
                threads[t] = new Thread(() -> {
                    long local = 0;
                    while (true) {
                        int idx = nextIdx.getAndIncrement();
                        if (idx >= primeLimit)
                            break;
                        local = addMod(local, branchFromPrime(idx), mod2);
                    }
                    partials[tid] = local;
                });
                threads[t].start();
            }

            long total = 0;
            for (int t = 0; t < cores; ++t) {
                try {
                    threads[t].join();
                    total = addMod(total, partials[t], mod2);
                } catch (InterruptedException e) {
                    e.printStackTrace();
                }
            }

            return total;
        }
    }

    public static String solve() {
        long N = 33557799775533L;
        Solver850 solver = new Solver850(N, MOD);
        long tMod = solver.computeTMod();
        long ans = (tMod - (tMod & 1L)) / 2L;
        return Long.toString(ans % MOD);
    }

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