Buggy python submission (verdict: Time Limit Exceeded):
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

is_prime = [False, True, False, False, False, True] * (max_A//6+1)
del is_prime[max_A+1:]
is_prime[1:4] = False, True, True
for i in range(5, int(max_A**0.5)+1):
    if is_prime[i]:
        is_prime[i*i::i] = [False] * (max_A//i-i+1)
primes = [p for p, is_p in enumerate(is_prime) if is_p]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]
g = [v * v % mod for v in g]
for p in primes:
    for k, g_kp in enumerate(g[p::p], 1):
        g[k] -= g_kp

modinv = [0, 1]
for i in range(2, max_A+1):
    modinv.append(mod - mod//i * modinv[mod%i] % mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv)))
ans %= mod
ans -= sum(A)
ans *= modinv[2]
ans %= mod
print(ans)


Fixed submission:
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [False, True, False, False, False, True] * (n//6+1)
    del is_prime[n+1:]
    is_prime[1:4] = False, True, True
    for i in range(5, int(n**0.5)+1):
        if is_prime[i]:
            is_prime[i*i::i] = [False] * (n//i-i+1)
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))

    # for k in range(1, max_A//p+1):
    #     g[k] -= g[k*p]  # 包除原理ヤバい

    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] -= g_kp  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)

Buggy python submission (verdict: Time Limit Exceeded):
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
#max_A = 10**6
N = int(input())
A = list(map(int, input().split()))
max_A = max(A)

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [False, True, False, False, False, True] * (n//6+1)
    del is_prime[n+1:]
    is_prime[1:4] = False, True, True
    for i in range(5, int(n**0.5)+1):
        if is_prime[i]:
            is_prime[i*i::i] = [False] * (n//i-i+1)
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))

    # for k in range(1, max_A//p+1):
    #     g[k] -= g[k*p]  # 包除原理ヤバい

    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] -= g_kp  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)


Fixed submission:
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [False, True, False, False, False, True] * (n//6+1)
    del is_prime[n+1:]
    is_prime[1:4] = False, True, True
    for i in range(5, int(n**0.5)+1):
        if is_prime[i]:
            is_prime[i*i::i] = [False] * (n//i-i+1)
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))

    # for k in range(1, max_A//p+1):
    #     g[k] -= g[k*p]  # 包除原理ヤバい

    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] -= g_kp  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)

Buggy python submission (verdict: Time Limit Exceeded):
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [False, True, False, False, False, True] * (n//6+1)
    del is_prime[n+1:]
    is_prime[1:4] = False, True, True
    for i in range(5, int(n**0.5)+1):
        if is_prime[i]:
            is_prime[i*i::i] = [False] * (n//i-i+1)
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))

    # for k in range(1, max_A//p+1):
    #     g[k] -= g[k*p]  # 包除原理ヤバい

    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] -= g_kp  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)


Fixed submission:
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [False, True, False, False, False, True] * (n//6+1)
    del is_prime[n+1:]
    is_prime[1:4] = False, True, True
    for i in range(5, int(n**0.5)+1):
        if is_prime[i]:
            is_prime[i*i::i] = [False] * (n//i-i+1)
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))

    # for k in range(1, max_A//p+1):
    #     g[k] -= g[k*p]  # 包除原理ヤバい

    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] -= g_kp  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)

Buggy python submission (verdict: Time Limit Exceeded):
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [False, True, False, False, False, True] * (n//6+1)
    del is_prime[n+2:]
    is_prime[1:4] = False, True, True
    for i in range(5, int(n**0.5)+1):
        if is_prime[i]:
            for j in range(i*i, n+1, i):
                is_prime[j] = False
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))

    # for k in range(1, max_A//p+1):
    #     g[k] -= g[k*p]  # 包除原理ヤバい

    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] -= g_kp  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)


Fixed submission:
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [False, True, False, False, False, True] * (n//6+1)
    del is_prime[n+1:]
    is_prime[1:4] = False, True, True
    for i in range(5, int(n**0.5)+1):
        if is_prime[i]:
            is_prime[i*i::i] = [False] * (n//i-i+1)
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))

    # for k in range(1, max_A//p+1):
    #     g[k] -= g[k*p]  # 包除原理ヤバい

    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] -= g_kp  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)

Buggy python submission (verdict: Time Limit Exceeded):
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = sorted(map(int, input().split()))  # ソートしたらメモリアクセスがいい感じにならないかな？

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [True] * (n+1)
    is_prime[0] = False
    is_prime[1] = False
    for i in range(2, int(n**0.5)+1):
        if is_prime[i]:
            for j in range(i*i, n+1, i):
                is_prime[j] = False
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))

    # for k in range(1, max_A//p+1):
    #     g[k] -= g[k*p]  # 包除原理ヤバい

    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] -= g_kp  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)


Fixed submission:
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [False, True, False, False, False, True] * (n//6+1)
    del is_prime[n+1:]
    is_prime[1:4] = False, True, True
    for i in range(5, int(n**0.5)+1):
        if is_prime[i]:
            is_prime[i*i::i] = [False] * (n//i-i+1)
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))

    # for k in range(1, max_A//p+1):
    #     g[k] -= g[k*p]  # 包除原理ヤバい

    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] -= g_kp  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)

Buggy python submission (verdict: Time Limit Exceeded):
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [True] * (n+1)
    is_prime[0] = False
    is_prime[1] = False
    for i in range(2, int(n**0.5)+1):
        if is_prime[i]:
            for j in range(i*i, n+1, i):
                is_prime[j] = False
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))

    # for k in range(1, max_A//p+1):
    #     g[k] -= g[k*p]  # 包除原理ヤバい

    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] -= g_kp  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)


Fixed submission:
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [False, True, False, False, False, True] * (n//6+1)
    del is_prime[n+1:]
    is_prime[1:4] = False, True, True
    for i in range(5, int(n**0.5)+1):
        if is_prime[i]:
            is_prime[i*i::i] = [False] * (n//i-i+1)
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))

    # for k in range(1, max_A//p+1):
    #     g[k] -= g[k*p]  # 包除原理ヤバい

    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] -= g_kp  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)

Buggy python submission (verdict: Time Limit Exceeded):
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [True] * (n+1)
    is_prime[0] = False
    is_prime[1] = False
    for i in range(2, int(n**0.5)+1):
        if is_prime[i]:
            for j in range(i*i, n+1, i):
                is_prime[j] = False
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))
    for k in range(1, max_A//p+1):
        g[k] -= g[k*p]  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0] * (n+1)
    modinv[1] = 1
    for i in range(2, n+1):
        modinv[i] = mod - mod//i * modinv[mod%i] % mod
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)


Fixed submission:
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [False, True, False, False, False, True] * (n//6+1)
    del is_prime[n+1:]
    is_prime[1:4] = False, True, True
    for i in range(5, int(n**0.5)+1):
        if is_prime[i]:
            is_prime[i*i::i] = [False] * (n//i-i+1)
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))

    # for k in range(1, max_A//p+1):
    #     g[k] -= g[k*p]  # 包除原理ヤバい

    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] -= g_kp  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)

Buggy python submission (verdict: Time Limit Exceeded):
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [True] * (n+1)
    is_prime[0] = False
    is_prime[1] = False
    for i in range(2, int(n**0.5)+1):
        if is_prime[i]:
            for j in range(i*i, n+1, i):
                is_prime[j] = False
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))
    for k in range(1, max_A//p+1):
        g[k] -= g[k*p]  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0] * (n+1)
    modinv[1] = 1
    for i in range(2, n+1):
        q, r = divmod(mod, i)
        modinv[i] = mod - q * modinv[r] % mod
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum(gg * minv % mod for gg, minv in zip(g, modinv_list))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)


Fixed submission:
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [False, True, False, False, False, True] * (n//6+1)
    del is_prime[n+1:]
    is_prime[1:4] = False, True, True
    for i in range(5, int(n**0.5)+1):
        if is_prime[i]:
            is_prime[i*i::i] = [False] * (n//i-i+1)
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))

    # for k in range(1, max_A//p+1):
    #     g[k] -= g[k*p]  # 包除原理ヤバい

    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] -= g_kp  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)

Buggy python submission (verdict: Time Limit Exceeded):
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [True] * (n+1)
    is_prime[0] = False
    is_prime[1] = False
    for i in range(2, int(n**0.5)+1):
        if is_prime[i]:
            for j in range(i*i, n+1, i):
                is_prime[j] = False
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))
    for k in range(1, max_A//p+1):
        g[k] -= g[k*p]  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0] * (n+1)
    modinv[1] = 1
    for i in range(2, n+1):
        modinv[i] = mod - mod // i * modinv[mod%i] % mod
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum(gg * minv % mod for gg, minv in zip(g, modinv_list))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)


Fixed submission:
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [False, True, False, False, False, True] * (n//6+1)
    del is_prime[n+1:]
    is_prime[1:4] = False, True, True
    for i in range(5, int(n**0.5)+1):
        if is_prime[i]:
            is_prime[i*i::i] = [False] * (n//i-i+1)
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))

    # for k in range(1, max_A//p+1):
    #     g[k] -= g[k*p]  # 包除原理ヤバい

    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] -= g_kp  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)

Buggy python submission (verdict: Time Limit Exceeded):
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

import sys
input = sys.stdin.readline
mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [True] * (n+1)
    is_prime[0] = False
    is_prime[1] = False
    for i in range(2, int(n**0.5)+1):
        if is_prime[i]:
            for j in range(i*i, n+1, i):
                is_prime[j] = False
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))
    for k in range(1, max_A//p+1):
        g[k] -= g[k*p]  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0] * (n+1)
    modinv[1] = 1
    for i in range(2, n+1):
        modinv[i] = mod - mod // i * modinv[mod%i] % mod
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum(gg * minv for gg, minv in zip(g, modinv_list))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)


Fixed submission:
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [False, True, False, False, False, True] * (n//6+1)
    del is_prime[n+1:]
    is_prime[1:4] = False, True, True
    for i in range(5, int(n**0.5)+1):
        if is_prime[i]:
            is_prime[i*i::i] = [False] * (n//i-i+1)
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))

    # for k in range(1, max_A//p+1):
    #     g[k] -= g[k*p]  # 包除原理ヤバい

    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] -= g_kp  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)

Buggy python submission (verdict: Time Limit Exceeded):
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [True] * (n+1)
    is_prime[0] = False
    is_prime[1] = False
    for i in range(2, int(n**0.5) + 1):
        if is_prime[i]:
            for j in range(i*i, n+1, i):
                is_prime[j] = False
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))
    for k in range(1, max_A//p+1):
        g[k] -= g[k*p]  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0] * (n+1)
    modinv[1] = 1
    for i in range(2, n+1):
        modinv[i] = mod - mod // i * modinv[mod%i] % mod
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum(gg * minv for gg, minv in zip(g, modinv_list))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)


Fixed submission:
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [False, True, False, False, False, True] * (n//6+1)
    del is_prime[n+1:]
    is_prime[1:4] = False, True, True
    for i in range(5, int(n**0.5)+1):
        if is_prime[i]:
            is_prime[i*i::i] = [False] * (n//i-i+1)
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))

    # for k in range(1, max_A//p+1):
    #     g[k] -= g[k*p]  # 包除原理ヤバい

    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] -= g_kp  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)

Buggy python submission (verdict: Time Limit Exceeded):
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [True] * (n+1)
    is_prime[0] = False
    is_prime[1] = False
    for i in range(2, int(n**0.5) + 1):
        if is_prime[i]:
            for j in range(i*i, n+1, i):
                is_prime[j] = False
    return is_prime

is_prime = make_prime_checker(max_A)

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in range(2, max_A+1):
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    if is_prime[p]:
        for k in range(max_A//p, 0, -1):
            g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in range(2, max_A+1):
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))
    if is_prime[p]:
        for k in range(1, max_A//p+1):
            g[k] -= g[k*p]  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0] * (n+1)
    modinv[1] = 1
    for i in range(2, n+1):
        modinv[i] = mod - mod // i * modinv[mod%i] % mod
    return modinv

modinv_list = make_modinv_list(max_A)
ans = sum(gg * minv % mod for gg, minv in zip(g, modinv_list))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# ans = Σ_{i<j} lcm(a,b)

print(ans)


Fixed submission:
# 解説 AC

# 参考1: https://twitter.com/rickytheta/status/1175412019006074880
# 参考2: https://maspypy.com/atcoder-%E5%8F%82%E5%8A%A0%E6%84%9F%E6%83%B3-2019-09-21agc-038

mod = 998244353
max_A = 10**6
N = int(input())
A = list(map(int, input().split()))

# nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
def make_prime_checker(n):
    is_prime = [False, True, False, False, False, True] * (n//6+1)
    del is_prime[n+1:]
    is_prime[1:4] = False, True, True
    for i in range(5, int(n**0.5)+1):
        if is_prime[i]:
            is_prime[i*i::i] = [False] * (n//i-i+1)
    return is_prime

primes = [p for p, is_prime in enumerate(make_prime_checker(max_A)) if is_prime]  # max_A 以下の素数のリスト

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    # 倍数集合の高速ゼータ変換みたいなやつ  O((max_A)loglog(max_A))
    # 参考: http://noshi91.hatenablog.com/entry/2018/12/27/121649
    # 大量に約数列挙みたいなことをするときはこれで高速化できる場合が多そう？（みんぷろ 2018 本戦 A - Uncommon など）
    for k in range(max_A//p, 0, -1):
        g[k] += g[k*p]

# この時点で g[d] = Σ_{d|a} a  (A の要素のうち d の倍数であるものの総和)

g = [v * v % mod for v in g]

# この時点で
# g[d] = (Σ_{d|a} a)(Σ_{d|b} b)
#      = Σ_{(d|a)∧(d|b)} ab
#      = Σ_{d | gcd(a,b)} ab

# この式変形天才すぎないか？（これを使った一連の操作を gcd 畳み込み というらしい？）

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))

    # for k in range(1, max_A//p+1):
    #     g[k] -= g[k*p]  # 包除原理ヤバい

    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] -= g_kp  # 包除原理ヤバい

# この時点で g[d] = Σ_{gcd(a,b)=d} ab

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す O(n)
    modinv = [0, 1]
    for i in range(2, n+1):
        modinv.append(mod - mod//i * modinv[mod%i] % mod)
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
ans %= mod

# この時点で
# ans = Σ_d ( Σ_{gcd(a,b)=d} ab/d )
#     = Σ_d ( Σ_{gcd(a,b)=d} lcm(a,b) )
#     = Σ_a Σ_b lcm(a,b)

ans -= sum(A)
ans *= modinv_list[2]
ans %= mod

# この時点で ans = Σ_{i<j} lcm(a,b)

print(ans)