from random import shuffle,randrange
rd=randrange

class primes():
    def __init__(self, n):
        self.prime_num = n
        self.min_prime = [-1] * (self.prime_num + 1)  # 2以上の自然数に対して最小の素因数を表す
        self.min_prime[0] = 0
        self.min_prime[1] = 1
        i = 2
        self.prime = []
        self.memo_prifac = {}
        while i <= self.prime_num:
            if self.min_prime[i] == -1:
                self.min_prime[i] = i
                self.prime.append(i)
            for j in self.prime:
                if i * j > self.prime_num or j > self.min_prime[i]: break
                self.min_prime[j * i] = j
            i += 1

    def prifac(self, n):
        # 素因数分解した結果を返す
        if n in self.memo_prifac:
            return self.memo_prifac[n]
        res = {}
        x = n
        while x > 1:
            p = self.min_prime[x]
            if p in res:
                res[p] += 1
            else:
                res[p] = 1
            x //= p

        # self.memo_prifac[n] = res  #場合によってはこの行を消すと高速化

        return res

    def divisors(self, n):
        # 約数列挙 メモした方がいいかも
        if n== 1: return [1]
        prf = self.prifac(n)
        keys = [key for key in prf]

        def divsearch(i):
            if i == len(keys) - 1:
                return [keys[i] ** j for j in range(prf[keys[i]] + 1)]
            else:
                res = []
                subres = divsearch(i + 1)
                p = keys[i]
                for j in range(prf[p] + 1):
                    for node in subres:
                        res.append(node * p ** j)
                return res

        return divsearch(0)

pri=primes(10**6+100)
u=[0]*(10**6+100)
u[1]=1
for x in range(2,10**6+10):
    u[x]=1
    d=pri.prifac(x)
    for p in d:
        if d[p]>=2:u[x]=0
        u[x]*=-1


from math import gcd
mod=998244353
def naive(n,p):
    ans=0
    for bit in range(1,2**n):
        x=[]
        for i in range(n):
            if (bit>>i)&1:x.append(p[i])
        k=len(x)
        flag=1
        for i in range(k-1):
            if gcd(x[i],x[i+1])==1:flag=0
        ans+=flag
    return ans%mod

def sol1(n,P):
    p=[0]+P[:]
    ma=max(p)
    dp=[0]*(ma+1)
    for i in range(1,n+1):
        res=0
        for j in range(1,ma+1):
            if gcd(p[i],j)==1:res+=dp[j]
        dp[p[i]]+=sum(dp)+1-res
        dp[p[i]]%=mod
    return sum(dp)%mod

def sol2(n,P):
    p=[0]+P[:]
    for i in range(1,n+1):
        res=1
        for q in pri.prifac(p[i]):res*=q
        p[i]=res
    ma=max(p)
    dp=[0]*(ma+1)
    g=[0]*(ma+1)
    sumdp=0
    for i in range(1,n+1):
        f1=0
        div=[]
        for m in pri.divisors(p[i]):div.append(m)
        for m in div:
            f1+=u[m]*g[m]
        delta=sumdp+1-f1
        dp[p[i]]+=delta
        dp[p[i]]%=mod
        sumdp+=delta
        sumdp%=mod
        for m in div:
            g[m]+=delta
            g[m]%=mod
    return sum(dp)%mod

n=int(input())
a=list(map(int,input().split()))
print(sol2(n,a))


cnt=0
while 0:
    cnt+=1
    print(cnt)
    n=randrange(1,19)
    p=[rd(1,100) for i in range(n)]
    ansn=naive(n,p)
    ans1=sol1(n,p)
    ans2=sol2(n,p)
    if ans1!=ansn:
        print(n)
        print(*p)
        exit()