素因子分治 + 状压

注意到对于 $\le 3000$ 的数来说，要么全都是 $\le \sqrt{3000}$ 的素因子，要么是存在一个 $\ge \sqrt{3000}$ 的素因子，那么我们对于每个大素数开一个 2^16 的桶， $dp[i][bit]$ 表示第 i 个素数，小素数状态为 bit 时的欧拉函数和。

考虑新加入 a[i] 产生了哪些新贡献： a[i] 没有大素数的时候很好办，有大素数的时候要讨论：对于之前没有该素数的欧拉函数和，需要计算上该素数的 (p - 1) / p，而之前有过的就不需要。此外，大家都要 * 新加入进来的小素数。

为了不统计错误，我们必须对 a 进行排序，按照： ① 有无大素数 ② 大素数的值相同大素数的必须放到一起处理！！！不然就会多乘一些 (p - 1) / p！！！具体看代码：

void solve() {
    int n;
    cin >> n;
    vector<int> a(n + 1);
    for (int i = 1; i <= n; i++) {
        cin >> a[i];
    }

    sieve(3e3);

    auto get = [&](int x) -> int {
        vector<int> res;
        while (minp[x]) {
            int c = minp[x];
            res.emplace_back(c);
            while (x % c == 0) x /= c;
        }
        if (res.empty()) return 0;
        return *max_element(res.begin(), res.end());
    };

    sort(a.begin() + 1, a.end(), [&](int x, int y) {
        return get(x) < get(y);
    });
    // 必须必须必须！！！必须这么排序！！！

    vector<int> id(3005);

    vector<Z> P(primes.size());
    for (int i = 0; i < primes.size(); i++) {
        P[i] = Z(primes[i] - 1) / primes[i];
    }

    int sq = sqrt(3e3);
    int idx = 0;
    for (auto t : primes) {
        id[t] = idx++;
    }

    vector<Z> f(1 << 16, 1);
    for (int i = 0; i < (1 << 16); i++) {
        for (int j = 0; j < 16; j++) {
            if (i >> j & 1) {
                f[i] = f[i] * P[j];
            }
        }
    }

    vector dp(primes.size(), vector<Z>(1 << 16));
    dp[0][0] = 1;

    vector<Z> pre(1 << 16);
    pre[0] = 1;

    for (int i = 1; i <= n; i++) {
        int t = 0, pos = 0;
        int c = a[i], x;
        while (minp[c]) {
            x = minp[c];
            if (x > sq) {
                assert(t == 0);
                t = id[x];
            } else {
                pos |= 1 << id[x];
            }
            while (c % x == 0) c /= x;
        }
        for (int j = (1 << 16) - 1; j >= 0; j--) {
            int bit = (j | pos) ^ j;
            if (t) {

Z ndp = dp[t][j];
dp[t][j | pos] += (pre[j] - ndp) * a[i] * f[bit] * P[t]
               + ndp * a[i] * f[bit];

pre[j | pos] += (pre[j] - ndp) * a[i] * f[bit] * P[t]
             + ndp * a[i] * f[bit];

// 核心公式！！！大家都要 *f[bit]，其中一些情况要乘 P[t]

            } else {
                pre[j | pos] += dp[t][j] * a[i] * f[bit];
                dp[t][j | pos] += dp[t][j] * a[i] * f[bit];
            }
        }
    }

    Z ans = 0;
    for (int i = 0; i < primes.size(); i++) {
        for (int j = 0; j < (1 << 16); j++) {
            ans += dp[i][j];
        }
    }
    cout << ans << "\n";
}

void solve() {
    int n;
    cin >> n;
    vector<int> a(n + 1);
    for (int i = 1; i <= n; i++) {
        cin >> a[i];
    }

    sieve(1e3);

    auto get = [&](int x) -> int {
        vector<int> res;
        while (minp[x]) {
            int c = minp[x];
            res.emplace_back(c);
            while (x % c == 0) x /= c;
        }
        if (res.empty()) return 0;
        return *max_element(res.begin(), res.end());
    };

    sort(a.begin() + 1, a.end(), [&](int x, int y) {
        return get(x) < get(y);
    });

    vector<int> id(1005);

    int sq = sqrt(1e3);
    int idx = 0;
    for (auto t : primes) {
        id[t] = idx++;
    }

    int m = primes.size();
    vector dp(m, vector<Z>(1 << 11));
    dp[0][0] = 1;

    for (int i = 1; i <= n; i++) {
        int t = 0;
        int c = a[i], x;
        map<int, int> map;
        while (minp[c]) {
            x = minp[c];
            if (x > sq) {
                assert(t == 0);
                t = id[x];
            }
            while (c % x == 0) {
                if (x <= sq) map[id[x]]++;
                c /= x;
            }
        }
        auto ndp = dp[t];
        auto dp0 = dp[0];
        for (int j = (1 << 11) - 1; j >= 0; j--) {
            int bit = 0;
            Z add = 1;

            for (int k = 0; k < 11; k++) {
                int idx = k;
                int cnt = map[k];
                if (((j >> idx & 1) + cnt) & 1) {
                    bit |= 1 << idx;
                }
add *= power(Z(primes[idx]), (j >> idx & 1) + cnt >> 1);
            }

            if (t) {
                dp[t][bit] += dp0[j] * add;
                dp[0][bit] += ndp[j] * add * primes[t];
            } else {
                dp[0][bit] += dp0[j] * add;
            }
        }
    }
    cout << dp[0][0] - 1 << "\n";
}