description: Mr. Chanek has an array a of n integers. The prettiness value of a is denoted as:
  
  $$$∑_{i=1}^{n} {∑_{j=1}^{n} {\gcd(a_i, a_j) ⋅ \gcd(i, j)}}$$$
  
  where \gcd(x, y) denotes the greatest common divisor (GCD) of integers x and y.
  
  In other words, the prettiness value of an array a is the total sum of \gcd(a_i, a_j) ⋅ \gcd(i, j) for all pairs (i, j).
  
  Help Mr. Chanek find the prettiness value of a, and output the result modulo 10^9 + 7!
  
  Input
  
  The first line contains an integer n (2 ≤ n ≤ 10^5).
  
  The second line contains n integers a_1, a_2, …, a_n (1 ≤ a_i ≤ 10^5).
  
  Output
  
  Output an integer denoting the prettiness value of a modulo 10^9 + 7.
  
  Example
  
  Input
  
  
  5
  3 6 2 1 4
  
  
  Output
  
  
  77
solution: #include <bits/stdc++.h>
  using namespace std;
  using ll = long long;
  using pii = pair<int, int>;
  using pil = pair<int, ll>;
  using pli = pair<ll, int>;
  using pll = pair<ll, ll>;
  template <typename T>
  bool chmax(T &x, const T &y) {
    return (x < y) ? (x = y, true) : false;
  }
  template <typename T>
  bool chmin(T &x, const T &y) {
    return (x > y) ? (x = y, true) : false;
  }
  template <typename T>
  int flg(T x, int i) {
    return (x >> i) & 1;
  }
  template <typename T>
  void print(const vector<T> &v, T x = 0) {
    int n = v.size();
    for (int i = 0; i < n; i++) cout << v[i] + x << (i == n - 1 ? '\n' : ' ');
  }
  template <typename T>
  void printn(const vector<T> &v, T x = 0) {
    int n = v.size();
    for (int i = 0; i < n; i++) cout << v[i] + x << '\n';
  }
  template <typename T>
  int lb(const vector<T> &v, T x) {
    return lower_bound(begin(v), end(v), x) - begin(v);
  }
  template <typename T>
  int ub(const vector<T> &v, T x) {
    return upper_bound(begin(v), end(v), x) - begin(v);
  }
  template <typename T>
  void rearrange(vector<T> &v) {
    sort(begin(v), end(v));
    v.erase(unique(begin(v), end(v)), end(v));
  }
  template <typename T>
  vector<int> id_sort(const vector<T> &v, bool greater = false) {
    int n = v.size();
    vector<int> ret(n);
    iota(begin(ret), end(ret), 0);
    sort(begin(ret), end(ret),
         [&](int i, int j) { return greater ? v[i] > v[j] : v[i] < v[j]; });
    return ret;
  }
  struct io_setup {
    io_setup() {
      ios_base::sync_with_stdio(false);
      cin.tie(NULL);
      cout << fixed << setprecision(15);
    }
  } io_setup;
  const int inf = (1 << 30) - 1;
  const ll INF = (1LL << 60) - 1;
  const int MOD = 1000000007;
  template <int mod>
  struct Mod_Int {
    int x;
    Mod_Int() : x(0) {}
    Mod_Int(long long y) : x(y >= 0 ? y % mod : (mod - (-y) % mod) % mod) {}
    static int get_mod() { return mod; }
    Mod_Int &operator+=(const Mod_Int &p) {
      if ((x += p.x) >= mod) x -= mod;
      return *this;
    }
    Mod_Int &operator-=(const Mod_Int &p) {
      if ((x += mod - p.x) >= mod) x -= mod;
      return *this;
    }
    Mod_Int &operator*=(const Mod_Int &p) {
      x = (int)(1LL * x * p.x % mod);
      return *this;
    }
    Mod_Int &operator/=(const Mod_Int &p) {
      *this *= p.inverse();
      return *this;
    }
    Mod_Int &operator++() { return *this += Mod_Int(1); }
    Mod_Int operator++(int) {
      Mod_Int tmp = *this;
      ++*this;
      return tmp;
    }
    Mod_Int &operator--() { return *this -= Mod_Int(1); }
    Mod_Int operator--(int) {
      Mod_Int tmp = *this;
      --*this;
      return tmp;
    }
    Mod_Int operator-() const { return Mod_Int(-x); }
    Mod_Int operator+(const Mod_Int &p) const { return Mod_Int(*this) += p; }
    Mod_Int operator-(const Mod_Int &p) const { return Mod_Int(*this) -= p; }
    Mod_Int operator*(const Mod_Int &p) const { return Mod_Int(*this) *= p; }
    Mod_Int operator/(const Mod_Int &p) const { return Mod_Int(*this) /= p; }
    bool operator==(const Mod_Int &p) const { return x == p.x; }
    bool operator!=(const Mod_Int &p) const { return x != p.x; }
    Mod_Int inverse() const {
      assert(*this != Mod_Int(0));
      return pow(mod - 2);
    }
    Mod_Int pow(long long k) const {
      Mod_Int now = *this, ret = 1;
      for (; k > 0; k >>= 1, now *= now) {
        if (k & 1) ret *= now;
      }
      return ret;
    }
    friend ostream &operator<<(ostream &os, const Mod_Int &p) {
      return os << p.x;
    }
    friend istream &operator>>(istream &is, Mod_Int &p) {
      long long a;
      is >> a;
      p = Mod_Int<mod>(a);
      return is;
    }
  };
  using mint = Mod_Int<MOD>;
  int main() {
    int N;
    cin >> N;
    vector<int> a(N);
    for (int i = 0; i < N; i++) cin >> a[i];
    int M = 100000;
    vector<vector<int>> ds(M + 1);
    vector<ll> c(M + 1);
    iota(c.begin(), c.end(), 0);
    for (int i = 1; i <= M; i++) {
      ds[i].emplace_back(i);
      for (int j = 2 * i; j <= M; j += i) {
        c[j] -= c[i];
        ds[j].emplace_back(i);
      }
    }
    mint ans = 0;
    vector<ll> cnt(M + 1, 0);
    vector<bool> used(M + 1, false);
    for (int i = 1; i <= N; i++) {
      mint tmp = 0;
      for (int j = i - 1; j < N; j += i) {
        for (auto &e : ds[a[j]]) {
          used[e] = true;
          cnt[e]++;
        }
      }
      for (int j = i - 1; j < N; j += i) {
        for (auto &e : ds[a[j]]) {
          if (used[e]) {
            used[e] = false;
            tmp += c[e] * cnt[e] * cnt[e];
          }
          cnt[e]--;
        }
      }
      ans += tmp * c[i];
    }
    cout << ans << '\n';
  }
  