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

Solution (Python):


def Eratosthenes(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]]


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






SIZE=1000001; MOD=998244353 

SIZE += 1
inv = [0]*SIZE  
inv[1] = 1
for i in range(2,SIZE):
    inv[i] = MOD - (MOD//i)*inv[MOD%i]%MOD



import sys

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)

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

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)

Solution (Python):


SIZE=1000001; MOD=998244353 

SIZE += 1
inv = [0]*SIZE  
inv[1] = 1
for i in range(2,SIZE):
    inv[i] = MOD - (MOD//i)*inv[MOD%i]%MOD


import sys

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






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






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