Problem 448: Average Least Common Multiple

Project Euler Problem 448 Solution

EulerSolve provides an optimized solution for Project Euler Problem 448, Average Least Common Multiple, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Define $$A(n)=\frac{1}{n}\sum_{i=1}^{n}\operatorname{lcm}(i,n),\qquad S(N)=\sum_{k=1}^{N}A(k).$$ The goal is to evaluate $S(99999999019)\bmod 999999017$. A direct computation of every least common multiple is far too slow, so the solution rewrites the average in terms of multiplicative functions and then evaluates the resulting summatory functions with floor-division blocks. Mathematical Approach Step 1: Rewrite One Average by Grouping Equal GCDs Fix $n$. For each $i\in\{1,\dots,n\}$, let $d=\gcd(i,n)$, write $n=dk$, and write $i=dj$. Then $\gcd(j,k)=1$, and $$\operatorname{lcm}(i,n)=\frac{in}{\gcd(i,n)}=\frac{djn}{d}=nj.$$ So the contribution depends only on the reduced residue $j$ modulo $k$. For a fixed divisor $k\mid n$, the admissible $j$ are exactly the integers $1\le j\le k$ with $\gcd(j,k)=1$. Therefore $$\sum_{i=1}^{n}\operatorname{lcm}(i,n)=n\sum_{k\mid n}\sum_{\substack{1\le j\le k\\ \gcd(j,k)=1}} j.$$ For $k>1$, reduced residues come in pairs $j$ and $k-j$, so their sum is $$\sum_{\substack{1\le j\le k\\ \gcd(j,k)=1}} j=\frac{k\varphi(k)}{2}.$$ The special case $k=1$ contributes $1$....

Detailed mathematical approach

Problem Summary

Define

$$A(n)=\frac{1}{n}\sum_{i=1}^{n}\operatorname{lcm}(i,n),\qquad S(N)=\sum_{k=1}^{N}A(k).$$

The goal is to evaluate $S(99999999019)\bmod 999999017$. A direct computation of every least common multiple is far too slow, so the solution rewrites the average in terms of multiplicative functions and then evaluates the resulting summatory functions with floor-division blocks.

Mathematical Approach

Step 1: Rewrite One Average by Grouping Equal GCDs

Fix $n$. For each $i\in\{1,\dots,n\}$, let $d=\gcd(i,n)$, write $n=dk$, and write $i=dj$. Then $\gcd(j,k)=1$, and

$$\operatorname{lcm}(i,n)=\frac{in}{\gcd(i,n)}=\frac{djn}{d}=nj.$$

So the contribution depends only on the reduced residue $j$ modulo $k$. For a fixed divisor $k\mid n$, the admissible $j$ are exactly the integers $1\le j\le k$ with $\gcd(j,k)=1$. Therefore

$$\sum_{i=1}^{n}\operatorname{lcm}(i,n)=n\sum_{k\mid n}\sum_{\substack{1\le j\le k\\ \gcd(j,k)=1}} j.$$

For $k>1$, reduced residues come in pairs $j$ and $k-j$, so their sum is

$$\sum_{\substack{1\le j\le k\\ \gcd(j,k)=1}} j=\frac{k\varphi(k)}{2}.$$

The special case $k=1$ contributes $1$. Combining both cases gives the classical identity

$$\sum_{i=1}^{n}\operatorname{lcm}(i,n)=\frac{n}{2}\left(1+\sum_{k\mid n}k\varphi(k)\right),$$

hence

$$\boxed{A(n)=\frac{1}{2}\left(1+\sum_{k\mid n}k\varphi(k)\right).}$$

Step 2: Turn the Outer Sum into a Divisor Sum

Summing the formula for $A(n)$ over $1\le n\le N$ yields

$$S(N)=\frac{1}{2}\left(N+\sum_{n=1}^{N}\sum_{d\mid n}d\varphi(d)\right).$$

Swap the order of summation: every divisor $d$ contributes once for each multiple of $d$ up to $N$, that is, $\left\lfloor N/d\right\rfloor$ times. Define

$$R(N)=\sum_{d=1}^{N}d\varphi(d)\left\lfloor\frac{N}{d}\right\rfloor.$$

Then

$$\boxed{S(N)=\frac{N+R(N)}{2}.}$$

This is the key reduction: the original lcm-average problem becomes a weighted divisor summatory problem.

Step 3: Evaluate $R(N)$ by Quotient Blocks

Introduce the weighted totient prefix sum

$$P(x)=\sum_{n\le x}n\varphi(n).$$

When $\left\lfloor N/d\right\rfloor$ is constant on an interval $[l,r]$, we can collapse that whole interval into one prefix-difference:

$$R(N)=\sum_{[l,r]}\left\lfloor\frac{N}{l}\right\rfloor\bigl(P(r)-P(l-1)\bigr).$$

The intervals are determined by the standard rule $v=\left\lfloor N/l\right\rfloor$, $r=\left\lfloor N/v\right\rfloor$. There are only $O(\sqrt{N})$ distinct quotient values, so block decomposition removes the need to scan every $d\le N$ individually.

Step 4: Express $P(x)$ with the Möbius Function

Use the identity

$$\varphi(n)=\sum_{d\mid n}\mu(d)\frac{n}{d}.$$

Multiplying by $n$ gives

$$n\varphi(n)=\sum_{d\mid n}d\mu(d)\left(\frac{n}{d}\right)^2.$$

Now sum over $n\le x$, write $n=dt$, and separate the variables:

$$P(x)=\sum_{d\le x}d\mu(d)\sum_{t\le x/d}t^2.$$

Let

$$U(m)=\sum_{t=1}^{m}t^2=\frac{m(m+1)(2m+1)}{6}.$$

Then

$$\boxed{P(x)=\sum_{d\le x}d\mu(d)\,U\!\left(\left\lfloor\frac{x}{d}\right\rfloor\right).}$$

So the only missing ingredient is a fast way to evaluate the prefix sum of $d\mu(d)$.

Step 5: Recurrence for the Weighted Möbius Prefix

Define

$$M(x)=\sum_{n\le x}n\mu(n).$$

Set $f(n)=n\mu(n)$. The Dirichlet-convolution identity $\operatorname{id}*\mu=\varphi$ implies another useful relation:

$$\sum_{d\mid m} d\,f\!\left(\frac{m}{d}\right)=m\sum_{e\mid m}\mu(e)= \begin{cases} 1,&m=1,\\ 0,&m>1. \end{cases}$$

Summing this over $1\le m\le x$ and exchanging the order of summation gives

$$\sum_{d\le x} d\,M\!\left(\left\lfloor\frac{x}{d}\right\rfloor\right)=1.$$

Separating the term $d=1$ yields the recursion

$$\boxed{M(x)=1-\sum_{d=2}^{x} d\,M\!\left(\left\lfloor\frac{x}{d}\right\rfloor\right).}$$

Again, $\left\lfloor x/d\right\rfloor$ is constant on quotient blocks, so this becomes

$$M(x)=1-\sum_{[l,r]\subseteq[2,x]}\left(\sum_{d=l}^{r}d\right)M\!\left(\left\lfloor\frac{x}{l}\right\rfloor\right).$$

The arithmetic progression sum $\sum_{d=l}^{r}d=\frac{(l+r)(r-l+1)}{2}$ is what makes each block computable in constant time.

Worked Example: $n=10$ and $N=10$

For $n=10$, the divisors are $1,2,5,10$, and

$$1\cdot\varphi(1)=1,\qquad 2\cdot\varphi(2)=2,\qquad 5\cdot\varphi(5)=20,\qquad 10\cdot\varphi(10)=40.$$

Therefore

$$A(10)=\frac{1}{2}(1+1+2+20+40)=32.$$

For the full prefix $N=10$, compute

$$R(10)=\sum_{d=1}^{10}d\varphi(d)\left\lfloor\frac{10}{d}\right\rfloor=274,$$

$$S(10)=\frac{10+274}{2}=142.$$

This agrees with a direct brute-force evaluation of the first ten averages.

How the Code Works

The C++, Python, and Java implementations all follow the same plan. They precompute $\mu(n)$, $\varphi(n)$, and the small prefix sums of $n\mu(n)$ and $n\varphi(n)$ up to a threshold $L\approx N^{2/3}$ using a linear sieve. For arguments above that threshold, they evaluate the two large prefix functions recursively, cache every large result, and always group terms by equal floor-division quotients. The final summation for $R(N)$ uses the same block structure. Because the modulus is prime, the divisions by $2$ and $6$ are carried out through modular inverses.

Complexity Analysis

The sieve up to $L\approx N^{2/3}$ costs $O(L)$ time and $O(L)$ memory. Each uncached large query visits only the distinct intervals on which a floor quotient is constant, rather than every integer one by one. With memoization, the total amount of large-query work stays on the same practical scale as the preprocessing. For the target size $N\approx 10^{11}$, this reduces the problem from impossible brute force to an $O(N^{2/3})$-scale method in time and memory.

Footnotes and References

Problem page: https://projecteuler.net/problem=448
Euler's totient function: Wikipedia — Euler's totient function
Möbius function: Wikipedia — Möbius function
Dirichlet convolution: Wikipedia — Dirichlet convolution
Floor-division block technique: cp-algorithms — divisor summatory techniques

Problem 448 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <numeric>
#include <unordered_map>
#include <vector>

using int64 = long long;

namespace {

constexpr int64 MOD = 999999017LL;
constexpr int64 DEFAULT_N = 99999999019LL;

inline int64 mod_norm(int64 x) {
    x %= MOD;
    if (x < 0) x += MOD;
    return x;
}

inline int64 mod_mul_norm(int64 a, int64 b) {
    return static_cast<int64>((static_cast<uint64_t>(a) * static_cast<uint64_t>(b)) % MOD);
}

inline int64 mod_add(int64 a, int64 b) {
    a += b;
    if (a >= MOD) a -= MOD;
    if (a < 0) a += MOD;
    return a;
}

inline int64 mod_sub(int64 a, int64 b) {
    a -= b;
    if (a < 0) a += MOD;
    return a;
}

inline int64 mod_mul(int64 a, int64 b) {
    return mod_mul_norm(mod_norm(a), mod_norm(b));
}

int64 egcd(int64 a, int64 b, int64& x, int64& y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    int64 x1, y1;
    int64 g = egcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - y1 * (a / b);
    return g;
}

int64 mod_inv(int64 a) {
    int64 x, y;
    int64 g = egcd(a, MOD, x, y);
    if (g != 1) return 0;
    return mod_norm(x);
}

constexpr int64 INV2 = (MOD + 1) / 2;
const int64 INV6 = mod_inv(6);

int64 sum_arith(int64 l, int64 r) {
    if (l > r) return 0;
    int64 cnt = (r - l + 1) % MOD;
    int64 s = mod_norm(l + r);
    return mod_mul_norm(mod_mul_norm(s, cnt), INV2);
}

int64 sum_sq(int64 n) {
    n = mod_norm(n);
    int64 a = n;
    int64 b = mod_norm(n + 1);
    int64 c = mod_norm(2 * n + 1);
    return mod_mul_norm(mod_mul_norm(mod_mul_norm(a, b), c), INV6);
}

struct Solver {
    int64 N;
    int LIM = 0;

    std::vector<int> primes;
    std::vector<int8_t> mu;
    std::vector<uint32_t> phi;
    std::vector<uint8_t> is_comp;
    std::vector<int32_t> prefH;
    std::vector<int32_t> prefG;

    std::unordered_map<int64, int32_t> memoH;
    std::unordered_map<int64, int32_t> memoG;

    explicit Solver(int64 n) : N(n) {
        long double x = std::pow(static_cast<long double>(N), 2.0L / 3.0L);
        LIM = static_cast<int>(x + 10);
        if (LIM < 100) LIM = 100;

        sieve();

        memoH.reserve(1 << 20);
        memoG.reserve(1 << 20);
        memoH.max_load_factor(0.7f);
        memoG.max_load_factor(0.7f);
    }

    void sieve() {
        mu.assign(LIM + 1, 0);
        phi.assign(LIM + 1, 0);
        is_comp.assign(LIM + 1, 0);
        prefH.assign(LIM + 1, 0);
        prefG.assign(LIM + 1, 0);

        primes.clear();
        primes.reserve(LIM / 10);

        mu[1] = 1;
        phi[1] = 1;

        for (int i = 2; i <= LIM; ++i) {
            if (!is_comp[i]) {
                primes.push_back(i);
                mu[i] = -1;
                phi[i] = static_cast<uint32_t>(i - 1);
            }
            for (int p : primes) {
                int64 v = static_cast<int64>(i) * p;
                if (v > LIM) break;
                is_comp[static_cast<size_t>(v)] = 1;
                if (i % p == 0) {
                    mu[static_cast<size_t>(v)] = 0;
                    phi[static_cast<size_t>(v)] = phi[i] * static_cast<uint32_t>(p);
                    break;
                }
                mu[static_cast<size_t>(v)] = static_cast<int8_t>(-mu[i]);
                phi[static_cast<size_t>(v)] = phi[i] * static_cast<uint32_t>(p - 1);
            }
        }

        for (int i = 1; i <= LIM; ++i) {
            int64 addH = 0;
            if (mu[i] == 1) {
                addH = i;
            } else if (mu[i] == -1) {
                addH = MOD - i;
            }
            prefH[i] = static_cast<int32_t>(mod_add(prefH[i - 1], addH));

            int64 addG = (static_cast<int64>(i) * static_cast<int64>(phi[i])) % MOD;
            prefG[i] = static_cast<int32_t>(mod_add(prefG[i - 1], addG));
        }
    }

    int32_t H(int64 n) {
        if (n <= 0) return 0;
        if (n <= LIM) return prefH[static_cast<size_t>(n)];
        auto it = memoH.find(n);
        if (it != memoH.end()) return it->second;

        int64 res = 1 % MOD;
        int64 l = 2;
        while (l <= n) {
            int64 v = n / l;
            int64 r = n / v;
            int64 coef = sum_arith(l, r);
            res = mod_sub(res, mod_mul_norm(coef, H(v)));
            l = r + 1;
        }

        int32_t out = static_cast<int32_t>(res);
        memoH.emplace(n, out);
        return out;
    }

    int32_t G(int64 n) {
        if (n <= 0) return 0;
        if (n <= LIM) return prefG[static_cast<size_t>(n)];
        auto it = memoG.find(n);
        if (it != memoG.end()) return it->second;

        int64 res = 0;
        int64 l = 1;
        while (l <= n) {
            int64 v = n / l;
            int64 r = n / v;
            int64 mu_seg = mod_sub(H(r), H(l - 1));
            res = mod_add(res, mod_mul_norm(mu_seg, sum_sq(v)));
            l = r + 1;
        }

        int32_t out = static_cast<int32_t>(res);
        memoG.emplace(n, out);
        return out;
    }

    int64 T(int64 n) {
        int64 res = 0;
        int64 l = 1;
        while (l <= n) {
            int64 v = n / l;
            int64 r = n / v;
            int64 seg = mod_sub(G(r), G(l - 1));
            res = mod_add(res, mod_mul_norm(v % MOD, seg));
            l = r + 1;
        }
        return res;
    }

    int64 S(int64 n) {
        int64 ans = mod_add(mod_norm(n), T(n));
        return mod_mul(ans, INV2);
    }

    static int64 lcm_ll(int64 a, int64 b) {
        return a / std::gcd(a, b) * b;
    }

    static int64 brute_S(int n) {
        int64 total = 0;
        for (int k = 1; k <= n; ++k) {
            int64 sum = 0;
            for (int i = 1; i <= k; ++i) sum += lcm_ll(k, i);
            if (sum % k != 0) {
                std::cerr << "[VALIDATION] lcm-sum not divisible by k=" << k << "\n";
                std::exit(1);
            }
            total += sum / k;
        }
        return total;
    }

    void run_validations() {
        int64 brute_100 = brute_S(100);
        if (brute_100 != 122726LL) {
            std::cerr << "[VALIDATION] brute S(100)=" << brute_100 << " expected 122726\n";
            std::exit(1);
        }

        const int checks[] = {1, 2, 3, 10, 50, 100, 200};
        for (int n : checks) {
            int64 brute = brute_S(n) % MOD;
            int64 fast = S(n);
            if (brute != fast) {
                std::cerr << "[VALIDATION] mismatch S(" << n << "): brute=" << brute
                          << " fast=" << fast << "\n";
                std::exit(1);
            }
        }

        int max_check = std::min(LIM, 200000);
        int64 direct = 0;
        for (int i = 1; i <= max_check; ++i) {
            direct = mod_add(direct, mod_mul(static_cast<int64>(mu[i]), static_cast<int64>(i)));
            if (H(i) != direct) {
                std::cerr << "[VALIDATION] H(" << i << ") mismatch\n";
                std::exit(1);
            }
        }

        for (int i = 1; i <= max_check; i += 97) {
            if (G(i) != prefG[i]) {
                std::cerr << "[VALIDATION] G(" << i << ") mismatch\n";
                std::exit(1);
            }
        }
    }
};

}  // namespace

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    const int64 N = DEFAULT_N;
    Solver solver(N);
    solver.run_validations();
    std::cout << solver.S(N) << "\n";
    return 0;
}

Python

def solve():
    MOD = 999999017
    N = 99999999019

    def mod_norm(x): return x % MOD
    def mod_add(a, b): return (a + b) % MOD
    def mod_sub(a, b): return (a - b) % MOD
    def mod_mul(a, b): return (a % MOD) * (b % MOD) % MOD

    def mod_inv(a):
        a = a % MOD
        b = MOD
        x0, x1 = 1, 0
        while b:
            q = a // b
            a, b = b, a - q * b
            x0, x1 = x1, x0 - q * x1
        return x0 % MOD

    INV2 = (MOD + 1) // 2
    INV6 = mod_inv(6)

    def sum_arith(l, r):
        if l > r: return 0
        cnt = (r - l + 1) % MOD
        s = (l + r) % MOD
        return s * cnt % MOD * INV2 % MOD

    def sum_sq(n):
        n = n % MOD
        return n * ((n+1) % MOD) % MOD * ((2*n+1) % MOD) % MOD * INV6 % MOD

    # Build sieve
    LIM = int(N ** (2/3)) + 10
    if LIM < 100: LIM = 100

    mu = [0] * (LIM + 1)
    phi = [0] * (LIM + 1)
    is_comp = bytearray(LIM + 1)
    primes = []
    mu[1] = 1; phi[1] = 1
    for i in range(2, LIM + 1):
        if not is_comp[i]:
            primes.append(i)
            mu[i] = -1
            phi[i] = i - 1
        for p in primes:
            v = i * p
            if v > LIM: break
            is_comp[v] = 1
            if i % p == 0:
                mu[v] = 0
                phi[v] = phi[i] * p
                break
            mu[v] = -mu[i]
            phi[v] = phi[i] * (p - 1)

    prefH = [0] * (LIM + 1)
    prefG = [0] * (LIM + 1)
    for i in range(1, LIM + 1):
        addH = 0
        if mu[i] == 1: addH = i
        elif mu[i] == -1: addH = MOD - i
        prefH[i] = mod_add(prefH[i-1], addH)
        addG = (i * phi[i]) % MOD
        prefG[i] = mod_add(prefG[i-1], addG)

    memoH = {}
    memoG = {}

    def H(n):
        if n <= 0: return 0
        if n <= LIM: return prefH[n]
        if n in memoH: return memoH[n]
        res = 1 % MOD
        l = 2
        while l <= n:
            v = n // l
            r = n // v
            coef = sum_arith(l, r)
            res = mod_sub(res, coef * H(v) % MOD)
            l = r + 1
        memoH[n] = res
        return res

    def G(n):
        if n <= 0: return 0
        if n <= LIM: return prefG[n]
        if n in memoG: return memoG[n]
        res = 0
        l = 1
        while l <= n:
            v = n // l
            r = n // v
            mu_seg = mod_sub(H(r), H(l-1))
            res = mod_add(res, mu_seg * sum_sq(v) % MOD)
            l = r + 1
        memoG[n] = res
        return res

    def T(n):
        res = 0
        l = 1
        while l <= n:
            v = n // l
            r = n // v
            seg = mod_sub(G(r), G(l-1))
            res = mod_add(res, (v % MOD) * seg % MOD)
            l = r + 1
        return res

    ans = mod_add(N % MOD, T(N))
    ans = ans * INV2 % MOD
    return str(ans)

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

Java

import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

public class Euler448 {
    static final long MOD = 999999017L;
    static final long DEFAULT_N = 99999999019L;

    static long modNorm(long x) {
        x %= MOD;
        if (x < 0)
            x += MOD;
        return x;
    }

    static long modAdd(long a, long b) {
        long res = a + b;
        if (res >= MOD)
            res -= MOD;
        if (res < 0)
            res += MOD;
        return res;
    }

    static long modSub(long a, long b) {
        long res = a - b;
        if (res < 0)
            res += MOD;
        return res;
    }

    static long egcd(long a, long b, long[] xy) {
        if (b == 0) {
            xy[0] = 1;
            xy[1] = 0;
            return a;
        }
        long[] xy1 = new long[2];
        long g = egcd(b, a % b, xy1);
        xy[0] = xy1[1];
        xy[1] = xy1[0] - xy1[1] * (a / b);
        return g;
    }

    static long modInv(long a) {
        long[] xy = new long[2];
        long g = egcd(a, MOD, xy);
        if (g != 1)
            return 0;
        return modNorm(xy[0]);
    }

    static final long INV2 = (MOD + 1) / 2;
    static final long INV6 = modInv(6);

    static long sumArith(long l, long r) {
        if (l > r)
            return 0;
        long cnt = (r - l + 1) % MOD;
        long s = modNorm(l + r);
        return ((s * cnt) % MOD * INV2) % MOD;
    }

    static long sumSq(long n) {
        n = modNorm(n);
        long a = n;
        long b = modNorm(n + 1);
        long c = modNorm(2 * n + 1);
        return ((((a * b) % MOD) * c) % MOD * INV6) % MOD;
    }

    static class Solver {
        long N;
        int LIM;

        int[] prefH;
        int[] prefG;

        Map<Long, Integer> memoH = new HashMap<>();
        Map<Long, Integer> memoG = new HashMap<>();

        Solver(long n) {
            this.N = n;
            double x = Math.pow(N, 2.0 / 3.0);
            LIM = (int) (x + 10);
            if (LIM < 100)
                LIM = 100;

            sieve();
        }

        void sieve() {
            byte[] mu = new byte[LIM + 1];
            int[] phi = new int[LIM + 1];
            boolean[] isComp = new boolean[LIM + 1];
            List<Integer> primes = new ArrayList<>(LIM / 10);

            mu[1] = 1;
            phi[1] = 1;

            for (int i = 2; i <= LIM; i++) {
                if (!isComp[i]) {
                    primes.add(i);
                    mu[i] = -1;
                    phi[i] = i - 1;
                }
                for (int p : primes) {
                    long v = (long) i * p;
                    if (v > LIM)
                        break;
                    isComp[(int) v] = true;
                    if (i % p == 0) {
                        mu[(int) v] = 0;
                        phi[(int) v] = phi[i] * p;
                        break;
                    }
                    mu[(int) v] = (byte) -mu[i];
                    phi[(int) v] = phi[i] * (p - 1);
                }
            }

            prefH = new int[LIM + 1];
            prefG = new int[LIM + 1];

            long hSum = 0;
            long gSum = 0;

            for (int i = 1; i <= LIM; i++) {
                long addH = 0;
                if (mu[i] == 1) {
                    addH = i;
                } else if (mu[i] == -1) {
                    addH = MOD - i;
                }
                hSum = modAdd(hSum, addH);
                prefH[i] = (int) hSum;

                long addG = ((long) i * phi[i]) % MOD;
                gSum = modAdd(gSum, addG);
                prefG[i] = (int) gSum;
            }
        }

        int H(long n) {
            if (n <= 0)
                return 0;
            if (n <= LIM)
                return prefH[(int) n];
            Integer cached = memoH.get(n);
            if (cached != null)
                return cached;

            long res = 1;
            long l = 2;
            while (l <= n) {
                long v = n / l;
                long r = n / v;
                long coef = sumArith(l, r);
                res = modSub(res, (coef * H(v)) % MOD);
                l = r + 1;
            }

            int out = (int) res;
            memoH.put(n, out);
            return out;
        }

        int G(long n) {
            if (n <= 0)
                return 0;
            if (n <= LIM)
                return prefG[(int) n];
            Integer cached = memoG.get(n);
            if (cached != null)
                return cached;

            long res = 0;
            long l = 1;
            while (l <= n) {
                long v = n / l;
                long r = n / v;
                long muSeg = modSub(H(r), H(l - 1));
                res = modAdd(res, (muSeg * sumSq(v)) % MOD);
                l = r + 1;
            }

            int out = (int) res;
            memoG.put(n, out);
            return out;
        }

        long T(long n) {
            long res = 0;
            long l = 1;
            while (l <= n) {
                long v = n / l;
                long r = n / v;
                long seg = modSub(G(r), G(l - 1));
                res = modAdd(res, ((v % MOD) * seg) % MOD);
                l = r + 1;
            }
            return res;
        }

        long S(long n) {
            long ans = modAdd(modNorm(n), T(n));
            return (ans * INV2) % MOD;
        }
    }

    public static String solve() {
        Solver solver = new Solver(DEFAULT_N);
        return Long.toString(solver.S(DEFAULT_N));
    }

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