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: Time Limit Exceeded

Solution code:
N=int(input())
A=list(map(int,input().split()))
mod=998244353

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

data={i:[0,0] for i in range(1,10**6+1)}
for i in range(N):
    data[A[i]][0]+=A[i]
    data[A[i]][1]+=A[i]**2

for i in range(1,10**6+1):
    for j in range(1,10**6//i+1):
        data[i][0]+=data[i*j][0]
        data[i][1]+=data[i*j][1]

ans={i:0 for i in data}
for i in range(10**6,0,-1):
    for j in range(1,10**6//i+1):
        ans[i]-=ans[i*j]
        ans[i]%=mod
    ans[i]+=(data[i][0]**2-data[i][1])//2
    ans[i]%=mod

res=0
for i in ans:
    res+=inverse[i]*ans[i]
    res%=mod

print(res)

Submission status: Time Limit Exceeded

Solution code:
from collections import Counter
MOD = 998244353
n = int(input())
a = list(map(int, input().split()))
c = Counter(a)

pws = [pow(x, 2, MOD) for x in range(1000001)]

res = [0 for _ in range(1000001)]
ans = 0

for i in range(1000000, 0, -1):
	sum_of_sq, sq_of_sum = 0, 0
	for j in range(i, 1000001, i):
		sum_of_sq += pws[j] * c[j]
		sum_of_sq %= MOD
		sq_of_sum += j * c[j]
		sq_of_sum %= MOD
	sq_of_sum = pow(sq_of_sum, 2, MOD)
	res[i] = (sq_of_sum - sum_of_sq) * pow(2*i, MOD-2, MOD)
	res[i] %= MOD
	for j in range(2*i, 1000001, i):
		res[i] -= res[j] * (j//i)
		res[i] %= MOD
	ans += res[i]
	ans %= MOD

print(ans)

Submission status: Time Limit Exceeded

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

Submission status: Time Limit Exceeded

Solution code:
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
l = [0] * C
for x in A:
    h[x] = (h[x] + x) % mod
    l[x] = (l[x] + x**2) % mod
for i in range(1, C):
    for j in range(i*2, C, i):
        h[i] = (h[i] + h[j]) % mod
        l[i] = (l[i] + l[j]) % mod
inv2 = inv_list[2]
g = [(h[i]**2 - l[i]) * inv2 % mod for i in range(C)]
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
print(ans)

Submission status: Time Limit Exceeded

Solution code:
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)
    
    #数列x[i]はi*(iの個数)
    for i in range(N):
        x[A[i]]+=A[i]
        
    #print(x)
    
    #上位集合のゼータ変換
    for k in range(1,MA+1):
        for s in range(k*2,MA+1,k):#sが2kからkずつ増える=>sはkの倍数(自分自身を数えないため，2kから)
            x[k]=(x[k]+x[s])%mod
    
    #print(x)
    #この段階でx[i]はiを約数に持つやつの和(個数も加味している)
    
    #積演算
    for i in range(MA+1):
        x[i]=(x[i]*x[i])%mod
        
    #print(x)
    
    
    #メビウス変換で戻す
    for k in range(MA,0,-1):
        for s in range(k*2,MA+1,k):
            x[k]=(x[k]-x[s])%mod
            
    #この段階でx[i]はiをgcdに持つもの同士をかけたものの和
    #print(x)
    
    #gcdでわる
    for i in range(2,MA+1):
        x[i]=(x[i]*pow(i,mod-2,mod))%mod
        
    #この段階でx[i]はiをgcdに持つもの同士のlcaの和
    #print(x)
    
    """
    for i in range(MA+1):
        if x[i]!=0:
            print(i,x[i])"""
    
    
    ans=0
    for i in range(1,MA+1):
        ans=(ans+x[i])%mod
        #xの和を取る
        
    for i in range(N):
        ans=(ans-A[i])%mod
        #A[i],A[i]の組を除去する．A[i]動詞のくみはlcaがA[i]になる
        
    ans=(ans*499122177)%mod#A[i]*A[j] と　A[j]*A[i]の被り除去
    
    print(ans)
        
    
    

main()

Submission status: Time Limit Exceeded

Solution code:
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)
    
    #数列x[i]はi*(iの個数)
    for i in range(N):
        x[A[i]]+=A[i]
        
    #print(x)
    
    #上位集合のゼータ変換
    for k in range(1,MA+1):
        for s in range(k*2,MA+1,k):#sが2kからkずつ増える=>sはkの倍数(自分自身を数えないため，2kから)
            x[k]=(x[k]+x[s])%mod
    
    #print(x)
    #この段階でx[i]はiを約数に持つやつの和(個数も加味している)
    
    #積演算
    for i in range(MA+1):
        x[i]=(x[i]*x[i])%mod
        
    #print(x)
    
    
    #メビウス変換で戻す
    for k in range(MA,0,-1):
        for s in range(k*2,MA+1,k):
            x[k]=(x[k]-x[s])%mod
            
    #この段階でx[i]はiをgcdに持つもの同士をかけたものの和
    #print(x)
    
    #gcdでわる
    for i in range(2,MA+1):
        x[i]=(x[i]*pow(i,mod-2,mod))%mod
        
    #この段階でx[i]はiをgcdに持つもの同士のlcaの和
    #print(x)
    
    """
    for i in range(MA+1):
        if x[i]!=0:
            print(i,x[i])"""
    
    
    ans=(sum(x))%mod
    #xの和を取るa
        
    ans=(ans-sum(A))%mod
    #A[i],A[i]の組を除去する．A[i]動詞のくみはlcaがA[i]になる
       
    """ 
    ans=(ans*pow(2,mod-2,mod))%mod#A[i]*A[j] と　A[j]*A[i]の被り除去
    print(pow(2,mod-2,mod))=499122177
    """
    ans=(ans*499122177)%mod
    
    print(ans)
        
    
    

main()

Submission status: Time Limit Exceeded

Solution code:
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)
    
    #数列x[i]はi*(iの個数)
    for i in range(N):
        x[A[i]]+=A[i]
        
    #print(x)
    
    #上位集合のゼータ変換
    for k in range(1,MA+1):
        for s in range(k*2,MA+1,k):#sが2kからkずつ増える=>sはkの倍数(自分自身を数えないため，2kから)
            x[k]=(x[k]+x[s])
    
    #print(x)
    #この段階でx[i]はiを約数に持つやつの和(個数も加味している)
    
    #積演算
    for i in range(MA+1):
        x[i]=(x[i]*x[i])%mod
        
    #print(x)
    
    
    #メビウス変換で戻す
    for k in range(MA,0,-1):
        for s in range(k*2,MA+1,k):
            x[k]=(x[k]-x[s])%mod
            
    #この段階でx[i]はiをgcdに持つもの同士をかけたものの和
    #print(x)
    
    #gcdでわる
    for i in range(2,MA+1):
        x[i]=(x[i]*pow(i,mod-2,mod))%mod
        
    #この段階でx[i]はiをgcdに持つもの同士のlcaの和
    #print(x)
    
    """
    for i in range(MA+1):
        if x[i]!=0:
            print(i,x[i])"""
    
    
    ans=0
    for i in range(1,MA+1):
        ans=(ans+x[i])%mod
        #xの和を取る
        
    for i in range(N):
        ans=(ans-A[i])%mod
        #A[i],A[i]の組を除去する．A[i]動詞のくみはlcaがA[i]になる
        
    ans=(ans*pow(2,mod-2,mod))%mod#A[i]*A[j] と　A[j]*A[i]の被り除去
    
    print(ans)
        
    
    

main()

Submission status: Time Limit Exceeded

Solution code:
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)
    
    #数列x[i]はi*(iの個数)
    for i in range(N):
        x[A[i]]+=A[i]
        
    #print(x)
    
    #上位集合のゼータ変換
    for k in range(1,MA+1):
        for s in range(k*2,MA+1,k):#sが2kからkずつ増える=>sはkの倍数(自分自身を数えないため，2kから)
            x[k]=(x[k]+x[s])%mod
    
    #print(x)
    #この段階でx[i]はiを約数に持つやつの和(個数も加味している)
    
    #積演算
    for i in range(MA+1):
        x[i]=(x[i]*x[i])%mod
        
    #print(x)
    
    
    #メビウス変換で戻す
    for k in range(MA,0,-1):
        for s in range(k*2,MA+1,k):
            x[k]=(x[k]-x[s])%mod
            
    #この段階でx[i]はiをgcdに持つもの同士をかけたものの和
    #print(x)
    
    #gcdでわる
    for i in range(2,MA+1):
        x[i]=(x[i]*pow(i,mod-2,mod))%mod
        
    #この段階でx[i]はiをgcdに持つもの同士のlcaの和
    #print(x)
    
    """
    for i in range(MA+1):
        if x[i]!=0:
            print(i,x[i])"""
    
    
    ans=0
    for i in range(1,MA+1):
        ans=(ans+x[i])%mod
        #xの和を取る
        
    for i in range(N):
        ans=(ans-A[i])%mod
        #A[i],A[i]の組を除去する．A[i]動詞のくみはlcaがA[i]になる
        
    ans=(ans*pow(2,mod-2,mod))%mod#A[i]*A[j] と　A[j]*A[i]の被り除去
    
    print(ans)
        
    
    

main()

Submission status: Time Limit Exceeded

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

Submission status: Time Limit Exceeded

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

Submission status: Time Limit Exceeded

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[:]
    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)

Submission status: Time Limit Exceeded

Solution code:
import sys
import numpy as np


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 になるような数列を作成
    dcm = np.array(solo_invs, dtype=np.int64)
    for i in range(1, len(dcm)):
        dcm[2 * i::i] -= dcm[i]
    dcm %= MOD
    return dcm


n, *aaa = map(int, sys.stdin.buffer.read().split())
MOD = 998244353
cnt = np.bincount(aaa)
LIMIT = len(cnt) - 1
facts, invs, solo_invs = prepare(LIMIT, MOD)
dcm = decompose_inverses(solo_invs, MOD)
idx = np.arange(LIMIT + 1, dtype=np.int64)
base1 = idx * cnt
base2 = idx ** 2 * cnt

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

Submission status: Time Limit Exceeded

Solution code:
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
    res[i]%=mod
    ans+=res[i]*pow(i, mod-2, mod)
    ans%=mod
print(ans*(pow(2, mod-2, mod))%mod)

Submission status: Time Limit Exceeded

Solution code:
import sys
from collections import Counter
from operator import mul
from functools import reduce
readline = sys.stdin.readline

#エラトステネスの篩と素因数分解

M = 10**6 + 3
prime = [0]*M
for i in range(2, M):
    if prime[i]:
        continue
    for j in range(i, M, i):
        if not prime[j]:
            prime[j] = i

def pd(x):
    C = Counter()
    while x > 1:
        C[prime[x]] += 1
        x //= prime[x]
    return C

MOD = 998244353
N = int(readline())
A = list(map(int, readline().split()))
C = Counter(A)
mA = max(A) + 1
table = [0]*mA
table2 = [0]*mA
T = [0]*mA
i2 = pow(2, MOD-2, MOD)

for k in A:
    pk = pd(k)
    lk = list(pk.keys())
    lv = [v+1 for v in pk.values()]
    ll = len(pk)
    for s in range(reduce(mul, lv, 1)):
        res = 1
        for idx in range(ll):
            res *= pow(lk[idx], s%lv[idx])
            s //= lv[idx]
        table[res] += k
        table2[res] += k*k

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)

Submission status: Time Limit Exceeded

Solution code:
from itertools import combinations
def compute_lcm(x, y):
   if x > y:
       greater = x
   else:
       greater = y
   while(True):
       if((greater % x == 0) and (greater % y == 0)):
           lcm = greater
           break
       greater += 1
   return lcm

a=int(input())
A=[int(x) for x in input().split()]
perm=combinations(A,2)
x=0
for i in perm:
    ele1=i[0]
    ele2=i[1]
    x+=compute_lcm(ele1,ele2)
print(x%998244353)

Submission status: Time Limit Exceeded

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

import fractions
 
def lcm(x, y):
    return (x * y) // fractions.gcd(x, y)

n = int(input())
a = list(map(int, input().split()))
 
lst = [
  lcm(a[i], a[j])%998244353 for i in range(n) for j in range(i + 1, n)
]
print(sum(lst)%998244353)

Submission status: Time Limit Exceeded

Solution code:
#C
#モジュール
import numpy as np
import fractions
#入力
N=int(input())
A=np.array(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]/fractions.gcd(A[i],A[j]))%998244353
        
print(x%998244353)

Submission status: Time Limit Exceeded

Solution code:
#C
#モジュール
import fractions
#入力
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]/fractions.gcd(A[i],A[j]))
        
print(x%998244353)

Submission status: Time Limit Exceeded

Solution code:
from fractions import gcd

def lcm(x,y):
  return x * y // gcd(x,y)

n=int(input())
a=list(map(int,input().split()))
ans = 0

for i in range(n-1):
  for j in range(i+1,n):
    ans = (ans + lcm(a[i],a[j])) % 998244353

print(ans)

Submission status: Time Limit Exceeded

Solution code:
import numpy as np


N=int(input())
A=np.array(list(map(int, input().split())))
Sum=0
def gcd(a,b):
	while b!=0:
		a,b=b,a%b
	return a
for i in range(N-1):
    for j in range(i+1,N):
        Sum+=A[i]*A[j]/gcd(A[i],A[j])%998244353
print(int(Sum%998244353))

Submission status: Time Limit Exceeded

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


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

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

    def sum_cross(arr):
        s = 0
        s2 = 0
        for a in arr:
            s += a
            s2 += a ** 2
            s %= mod
            s2 %= mod
        return ((((s ** 2) % mod - s2) % mod) * half) % mod

    ans = 0
    for i in range(1, maxA + 1):
        Ai = []
        for j in range(i, maxA + 1, i):
            for k in range(cnt[j]):
                Ai.append(j)
        ans += w[i] * sum_cross(Ai)
        ans %= mod
    print(ans)


if __name__ == '__main__':
    main()

Submission status: Time Limit Exceeded

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


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

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


    def sum_cross(arr):
        s = 0
        s2 = 0
        for a in arr:
            s += a
            s2 += a ** 2
            s %= mod
            s2 %= mod
        return ((((s ** 2) % mod - s2) % mod) * half) % mod


    ans = 0
    for i in range(1, maxA + 1):
        Ai = []
        for j in range(i, maxA + 1, i):
            for k in range(cnt[j]):
                Ai.append(j)
        ans += w[i] * sum_cross(Ai)
        ans %= mod

    print(ans)

if __name__ == '__main__':
    main()

Submission status: Time Limit Exceeded

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

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

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


def sum_cross(arr):
    s = 0
    s2 = 0
    for a in arr:
        s += a
        s2 += a ** 2
        s %= mod
        s2 %= mod
    return ((((s ** 2) % mod - s2) % mod) * half) % mod


ans = 0
for i in range(1, maxA + 1):
    Ai = []
    for j in range(i, maxA + 1, i):
        for k in range(cnt[j]):
            Ai.append(j)
    ans += w[i] * sum_cross(Ai)
    ans %= mod

print(ans)

Submission status: Time Limit Exceeded

Solution code:
#! /usr/bin/env python3

### @file lcm.py
### @brief
### @author Yusuke Matsunaga (松永 裕介)
###
### Copyright (C) 2019 Yusuke Matsunaga
### All rights reserved.


magic_number = 998244353

max_num = 1000001

n = int(input())

tmp_list = input().split()
A = [ int(tmp) for tmp in tmp_list ]

# 出現回数と最大値を求める．
h = [ 0 for i in range(max_num) ]
max_val = 0
for val in A :
    h[val] += 1
    if max_val < val :
        max_val = val

W = [ 1 for i in range(max_val + 1) ]
for i in range(1, max_val + 1) :
    for j in range(i * 2, max_val + 1, i) :
        W[j] -= W[i] * (j // i)

ans = 0
for x in range(1, max_val + 1) :
    sum = 0
    diff = 0
    for val1 in range(x, max_val + 1, x) :
        h1 = h[val1]
        v1 = h1 * val1
        sum += v1
        diff += v1 * val1
    val = (sum * sum - diff) // 2
    if val == 0 :
        continue

    ans += W[x] * val // x

ans = ans % magic_number

print(ans)

Submission status: Time Limit Exceeded

Solution code:
#! /usr/bin/env python3

### @file lcm.py
### @brief
### @author Yusuke Matsunaga (松永 裕介)
###
### Copyright (C) 2019 Yusuke Matsunaga
### All rights reserved.

from fractions import Fraction

magic_number = 998244353
max_num = 1000001
np_mark = [ False for i in range(max_num) ]
primes = []
for i in range(2, max_num) :
    if np_mark[i] :
        continue
    for j in range(i * 2, max_num, i) :
        np_mark[j] = True

for i in range(2, max_num) :
    if not np_mark[i] :
        primes.append(i)


def factoring(val) :
    ans = []
    for prime in primes :
        val1 = val
        c = 0
        while (val1 % prime) == 0 :
            c += 1
            val1 /= prime
        if c > 0 :
            ans.append( (prime, c) )
    return ans


n = int(input())

tmp_list = input().split()
A = [ int(tmp) for tmp in tmp_list ]

# 出現回数と最大値を求める．
h = dict()
max_val = 0
for val in A :
    if val not in h :
        h[val] = 0
    h[val] += 1
    if max_val < val :
        max_val = val

W = dict()
W[1] = Fraction(1, 1)

ans = Fraction(0)
for x in range(1, max_val + 1) :
    val = 0
    for val1 in range(x, max_val + 1, x) :
        if val1 not in h :
            continue
        h1 = h[val1]
        val2_sum = 0
        for val2 in range(val1 + x, max_val + 1, x) :
            if val2 in h :
                val2_sum += h[val2] * val2
        val += val1 * h1 * val2_sum
        if h1 >= 2 :
            val += val1 * val1 * ((h1 * (h1 - 1)) // 2)
    if val == 0 :
        continue

    # W[x] の計算を行う．
    #print('W[{}] = '.format(x), end = '')
    if x > 1 :
        fac = factoring(x)
        nf = len(fac)
        n_array = [ c for p, c in fac ]
        w = Fraction(1, x)
        c_array = [ 0 for i in range(nf) ]
        while True :
            index = 1
            for j in range(nf) :
                c = c_array[j]
                if c == 0 :
                    continue
                prime, dummy = fac[j]
                for k in range(c) :
                    index *= prime
            done = False
            for j in range(nf) :
                c_array[j] += 1
                if c_array[j] <= n_array[j] :
                    break
                c_array[j] = 0
            else :
                break
            w -= W[index]
        W[x] = w
    #print(W[x])
    ans += W[x] * Fraction(val)

ans = ans % magic_number

print(ans)

Submission status: Time Limit Exceeded

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()))

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)

Submission status: Time Limit Exceeded

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()))
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)

Submission status: Time Limit Exceeded

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 * 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: Time Limit Exceeded

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+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)