#include <bits/stdc++.h>
#define debug(x) cerr << #x << ": " << x << endl
#define debugArray(x, n)                           \
  for (long long hoge = 0; (hoge) < (n); ++(hoge)) \
  cerr << #x << "[" << hoge << "]: " << x[hoge] << endl
using namespace std;

template <int mod>
struct ModInt {
  int x;
  ModInt() : x(0) {}
  ModInt(int64_t y) : x(y >= 0 ? y % mod : (mod - (-y) % mod)) {}
  ModInt &operator+=(const ModInt &p) {
    if ((x += p.x) >= mod) x -= mod;
    return *this;
  }
  ModInt &operator-=(const ModInt &p) {
    if ((x += mod - p.x) >= mod) x -= mod;
    return *this;
  }
  ModInt &operator*=(const ModInt &p) {
    x = (int)(1LL * x * p.x % mod);
    return *this;
  }
  ModInt &operator/=(const ModInt &p) { return *this *= p.inverse(); }
  ModInt operator-() const { return ModInt() - *this; }
  ModInt operator+(const ModInt &p) const { return ModInt(*this) += p; }
  ModInt operator-(const ModInt &p) const { return ModInt(*this) -= p; }
  ModInt operator*(const ModInt &p) const { return ModInt(*this) *= p; }
  ModInt operator/(const ModInt &p) const { return ModInt(*this) /= p; }
  bool operator==(const ModInt &p) const { return x == p.x; }
  bool operator!=(const ModInt &p) const { return x != p.x; }
  ModInt inverse() const {
    int a = x, b = mod, u = 1, v = 0, t;
    while (b) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v);
    return ModInt(u);
  }
  ModInt pow(int64_t e) const {
    ModInt ret(1);
    for (ModInt b = *this; e; e >>= 1, b *= b)
      if (e & 1) ret *= b;
    return ret;
  }
  friend ostream &operator<<(ostream &os, const ModInt &p) { return os << p.x; }
  friend istream &operator>>(istream &is, ModInt &a) {
    int64_t t;
    is >> t;
    a = ModInt<mod>(t);
    return (is);
  }
  static int modulo() { return mod; }
};

namespace number_theory {
vector<int> primes;
void init(int n) {
  primes.push_back(2);
  bool is_prime[n + 1];
  fill(is_prime, is_prime + n + 1, true);
  for (long long p = 3; p <= n; p += 2)
    if (is_prime[p]) {
      primes.push_back(p);
      for (long long i = p * p; i <= n; i += p) is_prime[i] = false;
    }
}
// f -> g s.t. g(n) = sum_{m|n} f(m)
template <typename T>
void divisor_zeta(vector<T> &f) {
  int n = f.size();
  for (int p : primes) {
    if (p > n) break;
    for (int i = 1; p * i < n; i++) f[p * i] += f[i];
  }
}
// f -> h s.t. f(n) = sum_{m|n} h(m)
template <typename T>
void divisor_mobius(vector<T> &f) {
  int n = f.size();
  for (int p : primes) {
    if (p > n) break;
    for (int i = (n - 1) / p; i > 0; i--) f[p * i] -= f[i];
  }
}
// get table of Mobius function
template <typename T = int>
vector<T> get_mu(int n) {
  vector<T> f(n + 1, T(0));
  f[1] = T(1);
  divisor_mobius(f);
  return f;
}
// get table of Euler's totient function
template <typename T = int>
vector<T> get_phi(int n) {
  vector<T> f(n + 1);
  iota(f.begin(), f.end(), 0);
  divisor_mobius(f);
  return f;
}
// get table of number-of-divisors function
template <typename T = int>
vector<T> get_d(int n) {
  vector<T> f(n + 1, 1);
  divisor_zeta(f);
  return f;
}
// get table of sum-of-divisors function
template <typename T = int>
vector<T> get_sigma(int n) {
  vector<T> f(n + 1);
  iota(f.begin(), f.end(), 0);
  divisor_zeta(f);
  return f;
}
template <typename T>
static vector<T> lcm_convolution(vector<T> a, vector<T> b) {
  int n = max(a.size(), b.size());
  a.resize(n);
  b.resize(n);
  divisor_zeta(a);
  divisor_zeta(b);
  for (int i = 0; i < n; i++) a[i] *= b[i];
  divisor_mobius(a);
  return a;
}
// f -> g s.t. g(n) = sum_{n|m} f(m)
template <typename T>
static void multiple_zeta(vector<T> &f) {
  int n = f.size();
  for (int p : primes) {
    if (p > n) break;
    for (int i = (n - 1) / p; i > 0; i--) f[i] += f[p * i];
  }
}
// f -> h s.t. f(n) = sum_{n|m} h(m)
template <typename T>
static void multiple_mobius(vector<T> &f) {
  int n = f.size();
  for (int p : primes) {
    if (p > n) break;
    for (int i = 1; p * i < n; i++) f[i] -= f[p * i];
  }
}
template <typename T>
static vector<T> gcd_convolution(vector<T> a, vector<T> b) {
  int n = max(a.size(), b.size());
  a.resize(n);
  b.resize(n);
  multiple_zeta(a);
  multiple_zeta(b);
  for (int i = 0; i < n; i++) a[i] *= b[i];
  multiple_mobius(a);
  return a;
}
}  // namespace number_theory

signed main() {
  cin.tie(0);
  ios::sync_with_stdio(false);
  using Mint = ModInt<998244353>;
  using namespace number_theory;
  const int MAX_A = 1e6;
  init(MAX_A);
  int N;
  cin >> N;
  vector<Mint> c(MAX_A + 1, 0);
  Mint sum = 0;
  for (int i = 0; i < N; i++) {
    int A;
    cin >> A;
    c[A] += A;
    sum += A;
  }
  auto tmp = gcd_convolution(c, c);
  Mint ans = 0;
  for (int d = 1; d <= MAX_A; d++) {
    ans += tmp[d] / d;
  }
  ans = (ans - sum) / 2;
  cout << ans << endl;
  return 0;
}


Measured memory_bytes: 16340