Let ~DP_i~ be the number of ways to cover all the points that are a multiple of ~i~ without covering any other point. It's easy to see our answer to this question is going to be ~DP_1~. Also let ~cnt_i~ be the number of chips with periodicity ~i~.

We will update our ~DP~ from last to first. Imagine we want to find the value of ~DP_x~. It's obvious that we can only use chips with a periodicity that is a multiple of ~x~. One way to cover all points is to cover all the point with only one chip, in which case the answer will be ~cnt_x~.

Another way is to use the chips with periodicities strictly larger than, and multiples of, ~x~. let the greatest common divisor (gcd) of the periodicity of all the chips we are using be ~d \times x~ where ~d > 1~ since if ~d = 1~ then we would have a contradiction. If we fix ~d~ we can say for this particular ~d~ the answer would be ~(DP_{d \times x})^d~, since all the points can be grouped into ~d~ different non-overlapping groups such that, for each group, there will be ~DP_{d \times x}~ different ways to cover them.

Using the fact above we will use the inclusion–exclusion principle to calculate ~DP_x~. For every prime number ~p~ we will add ~(DP_{p \times x})^p~ to the answer.

Let's show an example now: if we have ~6x~ we counted it with both ~2x~ and also ~3x~ so we should subtract it.

Let's define a square free number as a number that, in its decomposition into primes, each prime has a power less or equal to one. Now for each square free number we will add it to the answer if it has an odd number of primes in its decomposition and we will subtract it if it has an even number of primes in its decomposition. And hence we will update all our DP values.

Final order: ~O(N + L \times logL)~

C++ Solution:

#include <bits/stdc++.h>
using namespace std;

#define pb push_back
#define all(c) (c).begin(), (c).end()
#define rall(c) (c).rbegin(), (c).rend()
#define fr first
#define sc second

typedef long long ll;
typedef long double ld;
typedef pair <int, int> pii;

const int MAXN = 1e6 + 10, MAXM = 1e2 + 10, MOD = 1e9 + 7;

int l[MAXN], dp[MAXN], val[MAXN], cnt[MAXN];
vector <int> square_free;


inline int pw (int a, int b) {
  int res = 1;
  while (b) {
    if (b & 1) res = 1ll * res * a % MOD;
    a = 1ll * a * a % MOD;
    b >>= 1;
  }
  return res;
}

int main() {
  ios_base::sync_with_stdio(false), cin.tie(NULL), cout.tie(NULL);
  fill(val, val + MAXN, 1);
  for (int i = 2; i < MAXN; i++) {
    if (val[i] == i) square_free.pb(i);
    if (cnt[i] > 0) continue;

    for (int j = i; j < MAXN; j += i) {
      val[j] *= i;
      cnt[j]++;
    }

    square_free.pb(i);
  }

  int n;
  cin >> n;
  for (int i = 0; i < n; i++) {cin >> l[i]; dp[l[i]]++;}

  for (int i = MAXN - 1; i >= 0; i--) {
    for (auto j: square_free) {
      if (i * j >= MAXN) break;

      int x = pw(dp[j * i], j);
      if (cnt[j] % 2) dp[i] = (dp[i] + x) % MOD;
      else dp[i] = (dp[i] + MOD - x) % MOD;
    }
  }

  cout << dp[1] << endl;
}