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

Solution (Python):
mod=998244353
N=int(input())
A=[int(i) for i in input().split()]
M=10**6
B=[0 for i in range(M+1)]
for i in range(N):
    B[A[i]]+=1
C=[0 for i in range(M+1)]
for d in range(1,M+1):
    for dk in range(d,M+1,d):
        C[d]+=B[dk]*dk
isPrime=[1 for i in range(M+1)]
isPrime[0]=0;isPrime[1]=0
for i in range(M+1):
    if isPrime[i]==0:
        continue
    for j in range(2*i,M+1,i):
        isPrime[j]=0
Mebius=[1 for i in range(M+1)]
Mebius[0]=0
for i in range(2,M+1):
    if isPrime[i]==0:
        continue
    for j in range(i,M+1,i):
        Mebius[j]*=-1
for i in range(2,M+1):
    for j in range(i*i,M+1,i*i):
        Mebius[j]=0

f=[0 for i in range(M+1)]
for d in range(1,M+1):
    for n in range(d,M+1,d):
        f[n]+=d*Mebius[d]
for n in range(1,M+1):
    f[n]=f[n]*pow(n,mod-2,mod)
ans=0
for d in range(1,M+1):
    ans+=f[d]*C[d]*C[d]
    ans%=mod
tmp=0
for i in range(N):
    tmp+=A[i]
submit=(ans-tmp)*pow(2,mod-2,mod)
print(submit%mod)

Solution (Python):
import sys
readline = sys.stdin.buffer.readline

def make_modinv_list(n, mod=10**9+7):
    inv_list = [0]*(n+1)
    inv_list[1] = 1
    for i in range(2, n+1):
        inv_list[i] = (mod - mod//i * inv_list[mod%i] % mod)
    return inv_list

mod = 998244353
N = int(readline())
A = list(map(int, readline().split()))
C = max(A) + 1
inv_list = make_modinv_list(C, mod)
h = [0] * C
for x in A:
    h[x] = (h[x] + x) % mod
for i in range(1, C):
    for j in range(i*2, C, i):
        h[i] = (h[i] + h[j]) % mod
g = [0] * C
for i in range(1, C):
    g[i] = h[i] * h[i] % mod
for i in range(C-1, 0, -1):
    for j in range(i*2, C, i):
        g[i] = (g[i] - g[j]) % mod
ans = 0
for i in range(1, C):
    ans = (ans + g[i] * inv_list[i]) % mod
ans -= sum(A)
print(ans * inv_list[2] % mod)

Solution (Python):
import sys
input = sys.stdin.readline

def I(): return int(input())
def MI(): return map(int, input().split())
def LI(): return list(map(int, input().split()))

def main():
    mod=998244353
    N=I()
    A=LI()
    MA=max(A)
    x=[0]*(MA+1)
    
    
    for i in range(N):
        x[A[i]]+=A[i]
        
    
    
    
    for k in range(1,MA+1):
        for s in range(k*2,MA+1,k):
            x[k]=(x[k]+x[s])%mod
    
    
    
    
    
    for i in range(MA+1):
        x[i]=(x[i]*x[i])%mod
        
    
    
    
    
    for k in range(MA,0,-1):
        for s in range(k*2,MA+1,k):
            x[k]=(x[k]-x[s])%mod
            
    
    
    
    
    ans=0
    for i in range(1,MA+1):
        if x[i]!=0:
            x[i]=(x[i]*pow(i,mod-2,mod))%mod
            ans=(ans+x[i])%mod
        
    
    
        
    ans=(ans-sum(A))%mod
    
       
    """ 
    ans=(ans*pow(2,mod-2,mod))%mod
    print(pow(2,mod-2,mod))=499122177
    """
    ans=(ans*499122177)%mod
    
    print(ans)
        
    
    

main()

Solution (Python):
import sys

readline = sys.stdin.readline

ns = lambda: readline().rstrip()
ni = lambda: int(readline().rstrip())
nm = lambda: map(int, readline().split())
nl = lambda: list(map(int, readline().split()))

def primes(n):
    is_prime = [True] * (n + 1)
    is_prime[0] = is_prime[1] = False
    for i in range(2, int((n+1)**0.5)+1):
        if is_prime[i]:
            for j in range(i *2, n + 1, i):
                is_prime[j] = False
    res = [i for i in range(n+1) if is_prime[i]]
    return res

def make_modinv_list(n, mod=10**9+7):
    inv_list = [0]*(n+1)
    inv_list[1] = 1
    for i in range(2, n+1):
        inv_list[i] = (mod - mod//i * inv_list[mod%i] % mod)
    return inv_list

def solve():
  mod = 998244353
  n = ni()
  a = nl()
  m = max(a)
  s = -sum(a) % mod
  l = [0]*(m+1)
  for x in a:
    l[x] += x
  a = make_modinv_list(m, mod)
  pr = primes(m) 
  for i in pr:
    for j in range(m//i, 0, -1):
      l[j] += l[j*i]
  for i in range(m+1):
    l[i] = l[i] * l[i] % mod
  for i in pr:
    for j in range(1, m//i + 1):
      l[j] = (l[j] - l[j*i]) % mod
  for i in range(1, m+1):
    if l[i]:
      s = (s + l[i] * a[i]) % mod
  print(s * a[2] % mod)
  return

solve()

Solution (Python):
import sys

readline = sys.stdin.readline

ns = lambda: readline().rstrip()
ni = lambda: int(readline().rstrip())
nm = lambda: map(int, readline().split())
nl = lambda: list(map(int, readline().split()))

def modinv(x, mod):
    a, b = x, mod
    u, v = 1, 0
    while b:
        t = a // b
        a -= t * b; a, b = b, a
        u -= t * v; u, v = v, u
    return u % mod

def solve():
  n = ni()
  a = nl()
  m = max(a)
  mod = 998244353
  l = [0]*(m+1)
  for x in a:
    l[x] += x
  sieve = [1]*(m+1)
  sieve[0] = sieve[1] = 0
  for i in range(2, m+1):
    if sieve[i]:
      for j in range(m//i, 0, -1):
        sieve[j*i] = 0
        l[j] += l[j*i]
      sieve[i] = 1
  for i in range(m+1):
    l[i] = l[i] * l[i] % mod
  for i in range(2, m+1):
    if sieve[i]:
      for j in range(1, m//i + 1):
        l[j] -= l[j*i]
  s = -sum(a)
  for i in range(1, m+1):
    if l[i]:
      s = (s + l[i] * modinv(i, mod)) % mod
  print(s * modinv(2, mod) % mod)
  return

solve()

Solution (Python):
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):
    
    n = len(solo_invs)
    dcm = solo_invs[:]
    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):
    mulsum = sum(count[d::d]) ** 2 - sum(double[d::d]) % MOD
    if mulsum:
        ans = (ans + dcm[d] * mulsum * inv2) % MOD
print(ans)

Solution (Python):
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):
    
    n = len(solo_invs)
    dcm = solo_invs[:]
    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)

Solution (Python):
import sys
input = sys.stdin.readline


n = int(input())
a = list(map(int, input().split()))
MOD = 998244353
max_a = max(a)

c = [0] * (max_a + 1)
for i in range(1, max_a + 1):
    c[i] = pow(i, MOD - 2, MOD)

for i in range(1, max_a):
    j = i
    while True:
        j += i
        if j <= max_a:
            c[j] -= c[i]
            c[j] %= MOD
        else:
            break

cnt = [0] * (max_a + 1)
for i in range(n):
   cnt[a[i]] += 1

ans = 0
for i in range(1, max_a + 1):
    tmp0 = 0
    tmp1 = 0
    j = i
    while True:
        if j <= max_a:
            tmp0 += cnt[j] * j
            tmp1 += cnt[j] * (j ** 2)
        else:
            break
        j += i
    tmp0 *= tmp0
    ans += (tmp0 - tmp1) * c[i]
    ans %= MOD

print((ans * pow(2, MOD - 2, MOD)) % MOD)

Solution (Python):
n = int(input())
a = list(map(int, input().split()))
MOD = 998244353

max_a = max(a)
c = [0] * (max_a + 1)
for i in range(1, max_a + 1):
    c[i] = pow(i, MOD - 2, MOD)

for i in range(1, max_a):
    j = i
    while True:
        j += i
        if j <= max_a:
            c[j] -= c[i]
            c[j] %= MOD
        else:
            break

cnt = [0] * (max_a + 1)
for i in range(n):
   cnt[a[i]] += 1

ans = 0
for i in range(1, max_a + 1):
    tmp0 = 0
    tmp1 = 0
    j = i
    while True:
        if j <= max_a:
            tmp0 += cnt[j] * j
            tmp1 += cnt[j] * (j ** 2)
        else:
            break
        j += i
    tmp0 *= tmp0
    ans += (tmp0 - tmp1) * c[i]
    ans %= MOD

print((ans * pow(2, MOD - 2, MOD)) % MOD)

Solution (Python):
mod = 998244353
n = int(input())
a = list(map(int, input().split()))
m = max(a)
cnt = [0]*(m+1)
for ai in a:
    cnt[ai]+=1
ans = 0
res = [0]*(m+1)
for i in range(1, m+1)[::-1]:
    sm = 0
    sq = 0
    for j in range(i, m+1, i):
        sm+=j*cnt[j] 
        sq+=cnt[j]*j*j
        res[i]-=res[j]
    res[i]+=sm*sm-sq
    ans+=res[i]//i
    ans%=mod
print(ans*(pow(2, mod-2, mod))%mod)

Solution (Python):
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 + 1
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

zeta_gcd(fs, primes)

for i, x in enumerate(fs):
    fs[i] = x**2

mobius_gcd(fs, primes)

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)

Solution (Python):
N=int(input())
A=list(map(int,input().split()))
V=max(A)
mod=998244353
data=[0]*(V+1)
data[1]=1
for i in range(2,V+1):
    data[i]=-(mod//i)*data[mod%i]%mod
for i in range(1,V+1):
    for j in range(2,V//i+1):
        data[i*j]=(data[i*j]-data[i])%mod
lst=[0]*(V+1)
for a in A:
    lst[a]+=1
mod_2=pow(2,mod-2,mod)
ans=0
for i in range(1,V+1):
    sum_1=0
    sum_2=0
    for j in range(1,V//i+1):
        zzz=i*j
        sum_1+=zzz*lst[zzz]
        sum_2+=zzz**2*lst[zzz]
    ans=(ans+(sum_1**2-sum_2)*data[i]*mod_2)%mod
print(ans)

Solution (Python):
import sys
readline = sys.stdin.readline    
MOD = 998244353
N = int(readline())
A = list(map(int, readline().split()))
mA = max(A) + 1
T = [0]*mA
i2 = pow(2, MOD-2, MOD)


table = [0]*mA
table2 = [0]*mA
for a in A:
    table[a] = (table[a] + a)%MOD
    table2[a] = (table2[a] + a*a)%MOD
prime = [True]*mA
for p in range(2, mA):
    if not prime[p]:
        continue
    for k in range((mA-1)//p, 0, -1):
        table[k] = (table[k] + table[k*p])%MOD
        table2[k] = (table2[k] + table2[k*p])%MOD
        prime[k*p] = False

T = [(x1*x1-x2)*i2%MOD for x1, x2 in zip(table, table2)]

ans = 0
coeff = [0]*mA
for i in range(1, mA):
    k = T[i]
    l = - coeff[i] + pow(i, MOD-2, MOD)
    ans = (ans+l*k)%MOD
    for j in range(2*i, mA, i):
        coeff[j] = (coeff[j] + l)%MOD
print(ans)

Solution (Python):
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 
        t=0 
        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))%MOD*mf.inv(2)%MOD
    print(ans%MOD)
resolve()

Solution (Python):
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 
        t=0 
        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()

Solution (Python):
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 
        t=0 
        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()

Solution (Python):
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)

Solution (Python):
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)

Solution (Python):
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()

Solution (Python):



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


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]

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    
    
    for k in range(max_A//p, 0, -1):
        g[k] = (g[k] + g[k*p]) % mod



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[k] - g_kp) % mod  






def make_modinv_list(n, mod=10**9+7):
    
    
    modinv = [0, 1]
    for a in range(2, n+1):
        modinv.append(- (mod//a) * modinv[mod%a] % 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 -= sum(A) 
ans *= modinv_list[2] 
ans %= mod



print(ans)

Solution (Python):
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)

Solution (Python):





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):
    
    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]  
    
    g = [0] * (max_A+1)
    for a in A:
        g[a] += a
    for p in primes:
        
        
        
        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_list = make_modinv_list(max_A, mod)
    ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
    ans %= mod
    
    
    
    
    
    
    ans -= sum(A)
    ans *= modinv_list[2]
    ans %= mod
    
    
    
    print(ans)

main()

Solution (Python):





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):
    
    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]  
    
    g = [0] * (max_A+1)
    for a in A:
        g[a] += a
    for p in primes:
        
        
        
        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_list = make_modinv_list(max_A, mod)
    ans = sum((gg * minv % mod for gg, minv in zip(g, modinv_list)))
    ans %= mod
    
    
    
    
    
    
    ans -= sum(A)
    ans *= modinv_list[2]
    ans %= mod
    
    
    
    print(ans)

main()

Solution (Python):
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]  

    g = [0] * (max_A+1)
    for a in A:
        g[a] += a
    for p in primes:
        
        
        
        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()

Solution (Python):





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


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]  

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    
    
    
    for k in range(max_A//p, 0, -1):
        g[k] = (g[k] + g[k*p]) % mod



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[k] - g_kp) % mod  



def make_modinv_list(n, mod=10**9+7):
    
    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 -= sum(A)
ans *= modinv_list[2]
ans %= mod



print(ans)

Solution (Python):





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


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]  

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    
    
    
    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  



def make_modinv_list(n, mod=10**9+7):
    
    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 -= sum(A)
ans *= modinv_list[2]
ans %= mod



print(ans)

Solution (Python):





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


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]  

g = [0] * (max_A+1)
for a in A:
    g[a] += a
for p in primes:
    
    
    
    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  



def make_modinv_list(n, mod=10**9+7):
    
    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 -= sum(A)
ans *= modinv_list[2]
ans %= mod



print(ans)

Solution (Python):
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)

Solution (Python):
P = 998244353
N = int(input())
A = [int(a) for a in input().split()]
V = max(A) + 1
X = [0] * V
for a in A:
    X[a] += 1

fa = [1]
for i in range(1, V):
    fa.append(fa[-1]*i%P)
fainv = [pow(fa[-1], P-2, P)]
for i in range(1, V)[::-1]:
    fainv.append(fainv[-1]*i%P)
fainv = fainv[::-1]
W = [1] + [fainv[i] * fa[i-1] % 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 = 0
    s2 = 0
    for i in range(k, V, k):
        if X[i]:
            s1 += i * X[i]
            s2 += i * i * X[i]
    ans = (ans + (s1*s1-s2) * W[k]) % P
print(ans * fainv[2] % P)

Solution (Python):
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



            
ANS=SUMLIST[1]*SUMLIST[1]



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

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

Solution (Python):
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 = [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 = 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()