#include <bits/stdc++.h>
#define rep(i,x,y) for (int i = x; i <= y; i++)
#define repd(i,x,y) for (int i = x; i >= y; i--)
#define mid ((l + r) >> 1)
#define lch (rt << 1)
#define rch (rt << 1 | 1)

using namespace std;

using u32 = unsigned;
using i64 = long long;
using u64 = unsigned long long;
 
// TODO: Dynamic ModInt
 
template<typename T>
constexpr T power(T a, u64 b) {
    T res {1};
    for (; b != 0; b /= 2, a *= a) {
        if (b % 2 == 1) {
            res *= a;
        }
    }
    return res;
}
 
template<u32 P>
constexpr u32 mulMod(u32 a, u32 b) {
    return 1ULL * a * b % P;
}
 
template<u64 P>
constexpr u64 mulMod(u64 a, u64 b) {
    u64 res = a * b - u64(1.L * a * b / P - 0.5L) * P;
    res %= P;
    return res;
}
 
template<typename U, U P>
struct ModIntBase {
public:
    constexpr ModIntBase() : x {0} {}
    
    template<typename T>
    requires std::integral<T>
    constexpr ModIntBase(T x_) : x {norm(x_ % T {P})} {}
    
    constexpr static U norm(U x) {
        if ((x >> (8 * sizeof(U) - 1) & 1) == 1) {
            x += P;
        }
        if (x >= P) {
            x -= P;
        }
        return x;
    }
    
    constexpr U val() const {
        return x;
    }
    
    constexpr ModIntBase operator-() const {
        ModIntBase res;
        res.x = norm(P - x);
        return res;
    }
    
    constexpr ModIntBase inv() const {
        return power(*this, P - 2);
    }
    
    constexpr ModIntBase &operator*=(const ModIntBase &rhs) & {
        x = mulMod<P>(x, rhs.val());
        return *this;
    }
    
    constexpr ModIntBase &operator+=(const ModIntBase &rhs) & {
        x = norm(x + rhs.x);
        return *this;
    }
    
    constexpr ModIntBase &operator-=(const ModIntBase &rhs) & {
        x = norm(x - rhs.x);
        return *this;
    }
    
    constexpr ModIntBase &operator/=(const ModIntBase &rhs) & {
        return *this *= rhs.inv();
    }
    
    friend constexpr ModIntBase operator*(ModIntBase lhs, const ModIntBase &rhs) {
        lhs *= rhs;
        return lhs;
    }
    
    friend constexpr ModIntBase operator+(ModIntBase lhs, const ModIntBase &rhs) {
        lhs += rhs;
        return lhs;
    }
    
    friend constexpr ModIntBase operator-(ModIntBase lhs, const ModIntBase &rhs) {
        lhs -= rhs;
        return lhs;
    }
    
    friend constexpr ModIntBase operator/(ModIntBase lhs, const ModIntBase &rhs) {
        lhs /= rhs;
        return lhs;
    }

    friend constexpr std::ostream &operator<<(std::ostream &os, const ModIntBase &a) {
        return os << a.val();
    }
    
    friend constexpr bool operator==(ModIntBase lhs, ModIntBase rhs) {
        return lhs.val() == rhs.val();
    }
    
    friend constexpr bool operator!=(ModIntBase lhs, ModIntBase rhs) {
        return lhs.val() != rhs.val();
    }
    
    friend constexpr bool operator<(ModIntBase lhs, ModIntBase rhs) {
        return lhs.val() < rhs.val();
    }
    
private:
    U x;
};
 
template<u32 P>
using ModInt = ModIntBase<u32, P>;
 
template<u64 P>
using ModInt64 = ModIntBase<u64, P>;
 
constexpr u32 P = 998244353;
using Z = ModInt<P>;

const vector<int> primes = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53};

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0);

    int n;
    cin >> n;

    vector<int> a(n);
    rep(i,0,n-1) {
        cin >> a[i];
    }

    auto getp = [&](int i) {
        for (auto p : primes) {
            while (i % p == 0) {
                i /= p;
            }
        }
        return i;
    };

    sort(a.begin(), a.end(), 
        [&](int i, int j) {
            return getp(i) < getp(j);
        });

    vector<int> d(n), s(n);
    rep(i,0,n-1) {
        d[i] = a[i];
        for (int j = 0; j < primes.size(); j++) {
            int p = primes[j];
            while (d[i] % p == 0) {
                d[i] /= p;
                s[i] |= 1 << j;
            }
        }
    }

    const int S = primes.size();
    vector<Z> mulp(1 << S, 1), mulp_(1 << S, 1);
    for (int i = 1; i < (1 << S); i++) {
        for (int j = 0; j < S; j++) {
            if (i >> j & 1) {
                mulp[i] *= primes[j];
                mulp_[i] *= primes[j] - 1;
            }
        }
    }

    vector<array<Z, 2>> f(1 << S);
    f[0][0] = 1;
    rep(i,0,n-1) {
        
        for (int j = (1 << S) - 1; j >= 0; j--) {
            int x = j & s[i];
            int y = s[i] ^ x;

            if (d[i] == 1) {
                f[j | s[i]][0] += f[j][0] * mulp[x] * mulp_[y];
            } else {
                f[j | s[i]][1] += (f[j][0] * (d[i] - 1) + f[j][1] * d[i]) * mulp[x] * mulp_[y];
            }
        }

        if (i < n - 1 && d[i + 1] != d[i]) {
            for (int j = 0; j < (1 << S); j++) {
                f[j][0] += f[j][1];
                f[j][1] = 0;
            }
        }
    }

    Z ans{0};
    for (int i = 0; i < (1 << S); i++) {
        ans += f[i][0] + f[i][1];
    }
    cout << ans << '\n';

    return 0;
}