Problem:


Score : 700 points

Problem Statement

We have an integer sequence of length N: A_0,A_1,\cdots,A_{N-1}.

Find the following sum (\mathrm{lcm}(a, b) denotes the least common multiple of a and b):

\sum_{i=0}^{N-2} \sum_{j=i+1}^{N-1} \mathrm{lcm}(A_i,A_j)

Since the answer may be enormous, compute it modulo 998244353.

Constraints

1 \leq N \leq 200000

1 \leq A_i \leq 1000000

All values in input are integers.

Input

Input is given from Standard Input in the following format:

N
A_0\ A_1\ \cdots\ A_{N-1}

Output

Print the sum modulo 998244353.

Sample Input 1

3
2 4 6

Sample Output 1

22

\mathrm{lcm}(2,4)+\mathrm{lcm}(2,6)+\mathrm{lcm}(4,6)=4+6+12=22.

Sample Input 2

8
1 2 3 4 6 8 12 12

Sample Output 2

313

Sample Input 3

10
356822 296174 484500 710640 518322 888250 259161 609120 592348 713644

Sample Output 3

353891724



Sample input:
3
2 4 6


Sample output:
22

Submission status: Accepted

Solution code:
import sys
sys.setrecursionlimit(2147483647)
INF=float("inf")
MOD=998244353
input=lambda:sys.stdin.readline().rstrip()

class modfact(object):
    def __init__(self,n):
        fact=[1]*(n+1); invfact=[1]*(n+1)
        for i in range(1,n+1): fact[i]=i*fact[i-1]%MOD
        invfact[n]=pow(fact[n],MOD-2,MOD)
        for i in range(n-1,-1,-1): invfact[i]=invfact[i+1]*(i+1)%MOD
        self.__fact=fact; self.__invfact=invfact

    def inv(self,n):
        assert(n>0)
        return self.__fact[n-1]*self.__invfact[n]%MOD

    def fact(self,n):
        return self.__fact[n]

    def invfact(self,n):
        return self.__invfact[n]

    def comb(self,n,k):
        if(k<0 or n<k): return 0
        return self.__fact[n]*self.__invfact[k]*self.__invfact[n-k]%MOD

    def perm(self,n,k):
        if(k<0 or n<k): return 0
        self.__fact[n]*self.__invfact[k]%MOD

def prime(n):
    if n<=1: return []
    S=[1]*(n+1)
    S[0]=0; S[1]=0
    for i in range(2,n):
        if(S[i]==0): continue
        for j in range(2*i,n+1,i): S[j]=0
    return [p for p in range(n+1) if(S[p])]

def resolve():
    n=int(input())
    A=list(map(int,input().split()))
    V=max(A)
    C=[0]*(V+1)
    for a in A: C[a]+=1

    P=prime(V)
    W=[1]*(V+1)
    for p in P:
        for i in range(p,V+1,p):
            W[i]*=(1-p)
    mf=modfact(V)
    for i in range(1,V+1):
        W[i]=(W[i]*mf.inv(i))%MOD

    ans=0
    for d in range(1,V+1):
        s=0 # 和(後に2乗する)
        t=0 # 2乗の和
        for i in range(d,V+1,d):
            s+=i*C[i]%MOD
            t+=(i**2)*C[i]%MOD
        ans+=W[d]*(s**2-t)*mf.inv(2)
        ans%=MOD
    print(ans)
resolve()

Submission status: Accepted

Solution code:
import sys
sys.setrecursionlimit(2147483647)
INF=float("inf")
MOD=998244353
input=lambda:sys.stdin.readline().rstrip()

class modfact(object):
    def __init__(self,n):
        fact=[1]*(n+1); invfact=[1]*(n+1)
        for i in range(1,n+1): fact[i]=i*fact[i-1]%MOD
        invfact[n]=pow(fact[n],MOD-2,MOD)
        for i in range(n-1,-1,-1): invfact[i]=invfact[i+1]*(i+1)%MOD
        self.__fact=fact; self.__invfact=invfact

    def inv(self,n):
        assert(n>0)
        return self.__fact[n-1]*self.__invfact[n]%MOD

    def fact(self,n):
        return self.__fact[n]

    def invfact(self,n):
        return self.__invfact[n]

    def comb(self,n,k):
        if(k<0 or n<k): return 0
        return self.__fact[n]*self.__invfact[k]*self.__invfact[n-k]%MOD

    def perm(self,n,k):
        if(k<0 or n<k): return 0
        self.__fact[n]*self.__invfact[k]%MOD

def prime(n):
    if n<=1: return []
    S=[1]*(n+1)
    S[0]=0; S[1]=0
    for i in range(2,n):
        if(S[i]==0): continue
        for j in range(2*i,n+1,i): S[j]=0
    return [p for p in range(n+1) if(S[p])]

def resolve():
    n=int(input())
    A=list(map(int,input().split()))
    V=max(A)
    C=[0]*(V+1)
    for a in A: C[a]+=1

    P=prime(V)
    W=[1]*(V+1)
    for p in P:
        for i in range(p,V+1,p):
            W[i]*=(1-p)
    mf=modfact(V)
    for i in range(1,V+1):
        W[i]=(W[i]*mf.inv(i))%MOD

    ans=0
    for d in range(1,V+1):
        s=0 # 和(後に2乗する)
        t=0 # 2乗の和
        for i in range(d,V+1,d):
            s+=i*C[i]
            t+=(i**2)*C[i]
        ans+=W[d]*(s**2-t)//2
        ans%=MOD
    print(ans)
resolve()

Submission status: Accepted

Solution code:
N=int(input());A=list(map(int, input().split()));z=[0]*(max(A)+1);m=998244353;r=0
for a in A:z[a]+=a
for i in range(1,len(z)):
  for j in range(2*i,len(z),i):z[i]+=z[j]
for i in range(1,len(z)):z[i]*=z[i]
for i in range(len(z)-1,0,-1):
  for j in range(2*i,len(z),i):z[i]-=z[j]
for a in A:z[a]-=a**2
for i in range(1,len(z)):
  z[i]//=i*2;r=(r+z[i])%m
print(r)

Submission status: Accepted

Solution code:
import sys
input = sys.stdin.buffer.readline
N = int(input())
A = list(map(int, input().split()))
zeta = [0] * (max(A) + 1)
mod = 998244353
res = 0
for a in A:
  zeta[a] += a
for i in range(1, len(zeta)):
  for j in range(2 * i, len(zeta), i):
    zeta[i] += zeta[j]
for i in range(1, len(zeta)):
  zeta[i] *= zeta[i]
for i in range(len(zeta) - 1, 0, -1):
  for j in range(2 * i, len(zeta), i):
    zeta[i] -= zeta[j]
for a in A:
  zeta[a] -= a ** 2
for i in range(1, len(zeta)):
  zeta[i] //= 2
for i in range(1, len(zeta)):
  zeta[i] //= i
  res = (res + zeta[i]) % mod
print(res)

Submission status: Accepted

Solution code:
N = int(input())
A = list(map(int, input().split()))
mod = 998244353
half = pow(2, mod - 2, mod)
maxA = max(A) + 1


def main():
    cnt = [0] * maxA
    for Ai in A:
        cnt[Ai] += 1

    w = [pow(i, mod-2, mod) for i in range(maxA)]
    for i in range(1, maxA):
        for j in range(i * 2, maxA, i):
            w[j] = (w[j] - w[i]) % mod

    ans = 0
    for i in range(1, maxA):
        s = s2 = 0
        for j in range(i, maxA, i):
            s = (s + (j * cnt[j])) % mod
            s2 = (s2 + (j * j * cnt[j])) % mod
        ans += w[i] * ((((s ** 2) % mod - s2) % mod) * half) % mod
        ans %= mod
    print(ans)


if __name__ == '__main__':
    main()

Submission status: Accepted

Solution code:
# https://atcoder.jp/contests/agc038/submissions/7671675
# を参考にした、というよりもはやそのもの

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

# max_A 以下の素数のリスト
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))
    # 大きい素数から再帰的に計算することで、約数列挙の手間を省く。賢すぎでは
    for k in range(max_A//p, 0, -1):
        g[k] = (g[k] + g[k*p]) % mod

# この時点で 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
# ただし重複込み

for p in primes:
    # 倍数集合の高速メビウス変換みたいなやつ  O((max_A)loglog(max_A))
    for k, g_kp in enumerate(g[p::p], 1):  # 高速化（ゼータ変換の方は途中で g が変化するので高速化できない）
        g[k] = (g[k] - g_kp) % mod  # 包除原理ヤバい
# 包除原理についてはg[d] = (dの倍数の個数)で考えて見ると何をしているのか分かりやすい
# 包除原理を適用可能な関数をモジュラ関数という
# 参考：http://www.misojiro.t.u-tokyo.ac.jp/~murota/lect-kisosuri/modularfunc081022.pdf

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

def make_modinv_list(n, mod=10**9+7):
    # 0 から n までの mod 逆元のリストを返す、有限体より逆元が存在 O(n)
    # 参考：http://wakabame.hatenablog.com/entry/2019/04/08/023634
    modinv = [0, 1]
    for a in range(2, n+1):
        modinv.append(- (mod//a) * modinv[mod%a] % mod) # 1/a = -q * 1/r
    return modinv

modinv_list = make_modinv_list(max_A, mod)
ans = sum((gg % mod * 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) # a = bのケースを抜く
ans *= modinv_list[2] # *1/2に相当
ans %= mod

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

print(ans)

Submission status: Accepted

Solution code:
import sys
readline = sys.stdin.readline
readlines = sys.stdin.readlines
sys.setrecursionlimit(10 ** 7)

import numpy as np

MOD = 998244353

N = int(readline())
A = np.array(readline().split(),np.int32)

U = 10**6+100
counter = np.bincount(A,minlength=U)

def cumprod(arr):
    L = len(arr); Lsq = int(L**.5+1)
    arr = np.resize(arr,Lsq**2).reshape(Lsq,Lsq)
    for n in range(1,Lsq):
        arr[:,n] *= arr[:,n-1]; arr[:,n] %= MOD
    for n in range(1,Lsq):
        arr[n] *= arr[n-1,-1]; arr[n] %= MOD
    return arr.ravel()[:L]

x = np.arange(U,dtype=np.int64); x[0] = 1
fact = cumprod(x)
x = np.arange(U,0,-1,dtype=np.int64); x[0] = pow(int(fact[-1]),MOD-2,MOD)
fact_inv = cumprod(x)[::-1]

inv = np.zeros(U,dtype=np.int64)
inv[1:] = fact[:-1] * fact_inv[1:] % MOD

"""
・f(n) = sum_{gcd(a,b) == 1}ab
・g(n) = sum_{n|gcd(a,b)}ab
"""

"""
・S(n) = sum_{n|a} a
・T(n) = sum_{n|a} a^2
"""

x = np.arange(U,dtype=np.int64)*counter%MOD
y = np.arange(U,dtype=np.int64)**2%MOD*counter%MOD
S = np.zeros(U,dtype=np.int64)
T = np.zeros(U,dtype=np.int64)
M = int(U**.5+1)+10
for d in range(1,M):
    S[d] = (x[d::d]).sum()
    T[d] = (y[d::d]).sum()
for d in range(1,M):
    S[M:(U-1)//d+1] += x[M*d:U:d]
    T[M:(U-1)//d+1] += y[M*d:U:d]
S%=MOD
T%=MOD

G = (S*S-T)%MOD*inv[2]%MOD

is_prime = np.zeros(U,dtype=np.int64)
is_prime[2] = 1
is_prime[3::2] = 1
for p in range(3,M,2):
    if p*p > U:
        break
    if is_prime[p]:
        is_prime[p*p::p+p] = 0
primes = np.where(is_prime)[0]
pi = is_prime.cumsum()

coef = np.ones(U,dtype=np.int64)
for p in primes[:pi[M]]:
    coef[p::p] *= (1-p)
for n in range(1,M):
    p = primes[pi[M]:pi[(U-1)//n]]
    coef[p*n] *= (1-p)

answer = (G*coef%MOD*inv%MOD).sum()%MOD
print(answer)

Submission status: Accepted

Solution code:
# 解説 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

def make_prime_checker(n):
    # nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
    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

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

def main():
    mod = 998244353
    max_A = 10**6
    N = int(input())
    A = list(map(int, input().split()))
    
    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
    
    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)

main()

Submission status: Accepted

Solution code:
# 解説 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

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

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

def main():
    mod = 998244353
    max_A = 10**6
    N = int(input())
    A = list(map(int, input().split()))
    
    # nまでの自然数が素数かどうかを表すリストを返す  O(nloglogn)
    
    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
    
    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)

main()

Submission status: Accepted

Solution code:
def main():
    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)
main()

Submission status: Accepted

Solution code:
# 解説 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] + g[k*p]) % mod

# この時点で 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[k] - g_kp) % mod  # 包除原理ヤバい

# この時点で 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 % mod * 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)

Submission status: Accepted

Solution code:
# 解説 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 % mod * 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)

Submission status: Accepted

Solution code:
P = 998244353
N = int(input())
A = [int(a) for a in input().split()]
V = 1000001
X = [0] * V
for a in A: X[a] += 1
W = [1] + [pow(i, P-2, P) for i in range(1, V)]
for k in range(1, V):
    for i in range(k*2, V, k):
        W[i] = (W[i] - W[k]) % P

ans = 0
for k in range(1, V):
    s1 = s2 = 0
    for i in range(k, V, k):
        s1 += i * X[i]
        s2 += i * i * X[i]
    ans = (ans + (s1*s1-s2) * W[k]) % P

print(ans * 499122177 % P)

Submission status: Accepted

Solution code:
import sys
input = sys.stdin.readline

N=int(input())
A=list(map(int,input().split()))
mod= 998244353
MAX=max(A)

COUNTLIST=[0]*(MAX+1)
for a in A:
    COUNTLIST[a]+=a

SUMLIST=[0]*(MAX+1)
CLIST=[0]*(MAX+1)

for i in range(1,MAX+1):
    score=0
    for j in range(i,MAX+1,i):
        score+=COUNTLIST[j]
    SUMLIST[i]=score

for i in range(1,MAX+1):
    CLIST[i]-=i-1
    ctemp=CLIST[i]
    for j in range(i*2,MAX+1,i):
        CLIST[j]-=j//i*ctemp

#print(SUMLIST)
#print(CLIST)
            
ANS=SUMLIST[1]*SUMLIST[1]

#print(ANS)

for i in range(2,MAX+1):
    ANS+=CLIST[i]*SUMLIST[i]*SUMLIST[i]*pow(i,mod-2,mod)
    #print(ANS%mod)

print((ANS%mod-SUMLIST[1])*pow(2,mod-2,mod)%mod)

Submission status: Accepted

Solution code:
def main():
    import sys
    input = sys.stdin.readline
    n = int(input())
    a = tuple(map(int,input().split()))
    mod = 998244353

    v = max(a)
    inv = [0]*(v+1)
    inv[1] = 1
    for i in range(2,v+1):
        inv[i] = mod - (mod//i)*inv[mod%i]%mod
    #w
    w = [1]*(v+1)
    for i in range(2,v+1):
        w[i] = (inv[i]-w[i])%mod
        for j in range(i*2,v+1,i):
            w[j] = (w[j] + w[i])%mod


    #res
    res = 0
    num = [0]*(v+1)
    for e in a:
        num[e] += 1


    for d in range(1,v+1):
        s = 0
        t = 0
        for j in range(d,v+1,d):
            s += num[j]*(j//d)
            t += (num[j]*(j//d)*(j//d))
        AA = ((((s**2-t)//2)%mod*d)%mod)*d%mod
        res = (res + w[d]*AA%mod)%mod

    print(res%mod)

if __name__ == '__main__':
    main()

Submission status: Accepted

Solution code:
def main():
    import sys
    input = sys.stdin.readline
    n = int(input())
    a = tuple(map(int,input().split()))
    mod = 998244353

    v = max(a)
    inv = [0]*(v+1)
    inv[1] = 1
    for i in range(2,v+1):
        inv[i] = mod - (mod//i)*inv[mod%i]%mod
    #w
    w = [1]*(v+1)
    for i in range(2,v+1):
        w[i] = (inv[i]-w[i])%mod
        for j in range(i*2,v+1,i):
            w[j] = (w[j] + w[i])%mod


    #res
    res = 0
    num = [0]*(v+1)
    for e in a:
        num[e] += 1


    for d in range(1,v+1):
        s = 0
        t = 0
        for j in range(d,v+1,d):
            s = (s + num[j]*(j//d))
            t = (t + (num[j]*(j//d))*(j//d))
        AA = ((((s**2-t)//2)%mod*d)%mod)*d%mod
        res = (res + w[d]*AA%mod)%mod

    print(res%mod)

if __name__ == '__main__':
    main()

Submission status: Accepted

Solution code:
# coding: utf-8
# Your code here!
def Eratosthenes(N): #N以下の素数のリストを返す
    N+=1
    is_prime_list = [True]*N
    m = int(N**0.5)+1
    for i in range(3,m,2):
        if is_prime_list[i]:
            is_prime_list[i*i::2*i]=[False]*((N-i*i-1)//(2*i)+1)
    return [2] + [i for i in range(3,N,2) if is_prime_list[i]]

#aを破壊的にgcd-zetaする
def zeta_gcd(a,primes):
    n = len(a)-1
    for p in primes:
        for i in range(n//p,0,-1):
            a[i] += a[p*i]
            a[i] %= MOD
#    return a

#aを破壊的にgcd-mobiusする
def mobius_gcd(a,primes):
    n = len(a)
    for p in primes:
        for i in range(1,n):
            if i*p >= n: break
            a[i] -= a[p*i] 
            a[i] %= MOD
#    return a

################################################

#################################################

SIZE=1000001; MOD=998244353 #ここを変更する

SIZE += 1
inv = [0]*SIZE  # inv[j] = j^{-1} mod MOD
inv[1] = 1
for i in range(2,SIZE):
    inv[i] = MOD - (MOD//i)*inv[MOD%i]%MOD

################################################
#################################################
import sys
#sys.setrecursionlimit(10**6)
readline = sys.stdin.readline 

n = int(input())
a = [int(i) for i in readline().split()]

num = [0]*SIZE
for i in a: num[i] += i

primes = Eratosthenes(SIZE)

zeta_gcd(num,primes)
ans = [i*i for i in num]
mobius_gcd(ans,primes)

res = 0
for i,c in enumerate(ans):
    res = (res + c*inv[i])%MOD
for i in a:
    res = (res-i)%MOD
        
print(res*inv[2]%MOD)

Submission status: Accepted

Solution code:
import sys
readline = sys.stdin.readline
readlines = sys.stdin.readlines
sys.setrecursionlimit(10 ** 7)

import numpy as np

MOD = 998244353

N = int(readline())
A = np.array(readline().split(),np.int32)

U = 10**6+100
counter = np.bincount(A,minlength=U)

def cumprod(arr):
    L = len(arr); Lsq = int(L**.5+1)
    arr = np.resize(arr,Lsq**2).reshape(Lsq,Lsq)
    for n in range(1,Lsq):
        arr[:,n] *= arr[:,n-1]; arr[:,n] %= MOD
    for n in range(1,Lsq):
        arr[n] *= arr[n-1,-1]; arr[n] %= MOD
    return arr.ravel()[:L]

x = np.arange(U,dtype=np.int64); x[0] = 1
fact = cumprod(x)
x = np.arange(U,0,-1,dtype=np.int64); x[0] = pow(int(fact[-1]),MOD-2,MOD)
fact_inv = cumprod(x)[::-1]

inv = np.zeros(U,dtype=np.int64)
inv[1:] = fact[:-1] * fact_inv[1:] % MOD

"""
・f(n) = sum_{gcd(a,b) == 1}ab
・g(n) = sum_{n|gcd(a,b)}ab
"""

"""
・S(n) = sum_{n|a} a
・T(n) = sum_{n|a} a^2
"""

x = np.arange(U,dtype=np.int64)*counter%MOD
y = np.arange(U,dtype=np.int64)**2%MOD*counter%MOD
S = np.zeros(U,dtype=np.int64)
T = np.zeros(U,dtype=np.int64)
M = int(U**.5+1)+10
for d in range(1,M):
    S[d] = (x[d::d]).sum()
    T[d] = (y[d::d]).sum()
for d in range(1,M):
    S[M:(U-1)//d+1] += x[M*d:U:d]
    T[M:(U-1)//d+1] += y[M*d:U:d]
S%=MOD
T%=MOD

G = (S*S-T)%MOD*inv[2]%MOD

coef = np.ones(U,dtype=np.int64)
for p in range(2,U):
    if coef[p] == 1:
        coef[p::p] *= (1-p)

answer = (G*coef%MOD*inv%MOD).sum()%MOD
print(answer)

Submission status: Accepted

Solution code:
# coding: utf-8
# Your code here!
SIZE=1000001; MOD=998244353 #ここを変更する

SIZE += 1
inv = [0]*SIZE  # inv[j] = j^{-1} mod MOD
inv[1] = 1
for i in range(2,SIZE):
    inv[i] = MOD - (MOD//i)*inv[MOD%i]%MOD


import sys
#sys.setrecursionlimit(10**6)
readline = sys.stdin.readline 

n = int(input())
a = [int(i) for i in readline().split()]

num = [0]*1000001
db = [0]*1000001
for i in a:
    num[i] += i
    db[i] += i*i
    db[i] %= MOD

def zeta(a,n):
    sieve = [True]*(n+1)
    for p in range(2,n):
        if sieve[p]:
            for k in range(n//p,0,-1):
                sieve[k*p] = False
                a[k] += a[k*p]
                a[k] %= MOD
    return a


def meb(a,n):
    sieve = [True]*(n+1)
    for p in range(2,n):
        if sieve[p]:
            for k in range(1,n):
                if k*p >= n: break
                sieve[k*p] = False
                a[k] -= a[k*p] 
                a[k] %= MOD
    return a


#n = 1000000

n = 1000000

#zn = [i*i for i in zeta(num,n)]
#zdb = zeta(db,n)
#print(zn[:20])
#print(zdb[:20])

zn = [i*i for i in zeta(num,n)]
zdb = zeta(db,n)
ans = [i-j for i,j in zip(zn, zdb)]

#ans = [i*i-j for i,j in zip(zeta(num,n), zeta(db,n))]

#print(ans[:20])
#print(meb(ans,n)[:20])

res = 0
for i,c in enumerate(meb(ans,n+1)):
    res = (res + c*inv[i])%MOD
        
print(res*inv[2]%MOD)

Submission status: Runtime Error

Solution code:
import sys


def prepare(n, MOD):
    f = 1
    factorials = [1]
    for m in range(1, n + 1):
        f *= m
        f %= MOD
        factorials.append(f)
    inv = pow(f, MOD - 2, MOD)
    invs = [1] * (n + 1)
    invs[n] = inv
    for m in range(n, 1, -1):
        inv *= m
        inv %= MOD
        invs[m - 1] = inv

    solo_invs = [0] + [f * i % MOD for f, i in zip(factorials, invs[1:])]

    return factorials, invs, solo_invs


def decompose_inverses(solo_invs, MOD):
    # 各整数 g に対して、g の約数である各 i について dcm[i] を全て足すと 1/g になるような数列を作成
    n = len(solo_invs)
    dcm = solo_invs.copy()
    for i in range(1, n):
        d = dcm[i]
        for j in range(2 * i, n, i):
            dcm[j] -= d
    for i in range(1, n):
        dcm[i] %= MOD
    return dcm


n, *aaa = map(int, sys.stdin.buffer.read().split())
MOD = 998244353
LIMIT = max(aaa)
count = [0] * (LIMIT + 1)
double = [0] * (LIMIT + 1)
for a in aaa:
    count[a] += a
    double[a] += a * a
_, _, solo_invs = prepare(LIMIT, MOD)
dcm = decompose_inverses(solo_invs, MOD)

ans = 0
inv2 = solo_invs[2]
for d in range(1, LIMIT + 1):
    ans = (ans + dcm[d] * (sum(count[d::d]) ** 2 - sum(double[d::d])) * inv2) % MOD
print(ans)

Submission status: Runtime Error

Solution code:
MOD = 998244353

def zeta_gcd(a,primes):
    n = len(a)-1
    for p in primes:
        for i in range(n//p,0,-1):
            a[i] += a[p*i]
            a[i] %= MOD
def mobius_gcd(a,primes):
    n = len(a)
    for p in primes:
        for i in range(1,n):
            if i*p >= n: break
            a[i] -= a[p*i] 
            a[i] %= MOD
def get_primes(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 not is_prime[i]:
            continue
        for j in range(i * 2, n + 1, i):
            is_prime[j] = False
    return [i for i in range(n + 1) if is_prime[i]]

NNN = 10**6
inverse = [0, 1]
for i in range( 2, NNN + 1 ):
    inverse.append( ( -inverse[MOD % i] * (MOD//i) ) % MOD )

N, = map(int, input().split())
X = list(map(int, input().split()))
primes = get_primes(NNN)
fs = [0]*NNN
for x in X:
    fs[x] += x
#print(fs)
zeta_gcd(fs, primes)
#print(fs)
for i, x in enumerate(fs):
    fs[i] = x**2
#print(fs)
mobius_gcd(fs, primes)
#print(fs)
r = 0
for i, f in enumerate(fs[1:]):
    i += 1
    r = (r + f*inverse[i]) % MOD
for x in X:
    r = (r - x) % MOD

print((r*inverse[2])%MOD)

Submission status: Runtime Error

Solution code:
N = int(input())
A = list(map(int,input().split()))
e = 998244353
import fractions
ans = 0
for i in range(N):
    x = A[i]
    for j in range(i+1,N):
        y = A[j]
        g = fractions.gcd(x,y)
        l = (x*y)/g
        s += l % e
print(s)

Submission status: Runtime Error

Solution code:
import math

n = int(input())
p = list(map(int, input().split()))

numsum = 0
for i in range(n-1):
    for j in  range(i+1,n):
        a = p[i]
        b = p[j]
        numsum+= a * b /math.gcd(a, b)

print(int(numsum%998244353))

Submission status: Runtime Error

Solution code:
from math import gcd
x=int(input())
a=list(map(int,input().split()))
lol=0
for i in range(x):
  for j in range(i+1,x):
    lol+=int(((a[i]*a[j])/gcd(a[i],a[j]))%998244353)
print(lol%998244353)

Submission status: Runtime Error

Solution code:
from math import gcd
x=int(input())
a=list(map(int,input().split()))
lol=0
for i in range(x):
  for j in range(i+1,x):
    lol+=((a[i]*a[j])/gcd(a[i],a[j]))%998244353
print(lol%998244353)

Submission status: Runtime Error

Solution code:
import sys
input = sys.stdin.buffer.readline
N = int(input())
A = list(map(int, input().split()))
zeta = [0] * (max(A) + 1)
mod = 998244353
for a in A:
  zeta[a] += a
for i in range(1, len(zeta)):
  for j in range(2 * i, len(zeta), i):
    zeta[i] += zeta[j]
for i in range(1, len(zeta)):
  zeta[i] *= zeta[i]
for i in range(len(zeta) - 1, 0, -1):
  for j in range(2 * i, len(zeta), i):
    zeta[i] -= zeta[j]
for a in A:
  zeta[a] -= a ** 2
for i in range(1, len(zeta)):
  zeta[i] //= 2
for i in range(1, len(zeta)):
  zeta[i] //= i
  res = (res + zeta[i]) % mod
print(res)

Submission status: Runtime Error

Solution code:
#C
#モジュール
import math
#入力
N=int(input())
A=list(map(int, input().split()))
x=0
#処理
for i in range(N-1):
    for j in range(i+1,N):
        x+=int(A[i]*A[j]/math.gcd(A[i],A[j]))
        
print(x%998244353)

Submission status: Runtime Error

Solution code:
#C
#モジュール
import math
#入力
N=int(input())
A=list(map(int, input().split()))
x=0
#処理
for i in range(N-1):
    for j in range(i+1,N):
        x+=A[i]*A[j]//math.gcd(A[i],A[j])
        
print(x%998244353)