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)



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)
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**0.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**0.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)

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**0.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**0.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)

+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() -for p in range(2, U): -    if coef[p] == 1: -        coef[p::p] *= 1 - p +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