Problem: 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.

Submission attempts (in order):
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





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

for n in range(U-1,0,-1):
    G[n] -= G[2*n:U:n].sum()
    G[n] %= MOD

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

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





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)

Verdicts: Time Limit Exceeded

Accepted