EZOI [NOIP模拟赛][莫比乌斯反演]

求$\sum_{i_1=1}^{a_1}$$\sum_{i_2=1}^{a_2}$$...$$\sum_{i_n=1}^{a_n}{gcd(i_1,i_2,...,i_n)}^{gcd(i_1,i_2,...,i_n)}$ (2≤n≤1000， 1≤ai≤100000)

先来考虑$n=2$的情况：

通过莫比乌斯反演可以得到：

$ans=$$\sum_{G=1}^{a_1}f(G)⌊\frac{a_1}{G}⌋⌊\frac{a_2}{G}⌋$ 不妨设$a_1≤a_2$

其中：$f(G)=\sum_{d|G}\mu(\frac{G}{d})d^d$

同理可以得出$n>2$的情况

$ans=$$\sum_{G=1}^{a_1}f(G)⌊\frac{a_1}{G}⌋⌊\frac{a_2}{G}⌋...⌊\frac{a_n}{G}⌋$ 不妨设$a_1≤a_{2...n}$

其中：$f(G)=\sum_{d|G}\mu(\frac{G}{d})d^d$

现在面临的问题就是要解决$f(G)$

由于$f(G)$函数中的求和是整除，不方便处理，不如换一个角度思考。枚举$d$，对$d$的所有倍数累加和，这样就能够很快地预处理出$f(G)$函数。

之后与普通的莫比乌斯反演相同

具体实现过程如下：

#include <cstdio>
#include <cstdlib>
#include <iostream>
using namespace std;

const int N = 100010, MOD = 1e9 + 7;
int a[N], mu[N], prime[N];
long long pre[N], f[N], po[N];
bool vis[N] = {1, 1};
int n, MAXN, pcnt, x, lim = 1 << 30, pos;
long long ans, res;

inline long long pow(long long a, long long b) {
    long long c = 1;
    while (b) {
        if (b & 1) {
            c = (ll)(c * a) % MOD; --b;
        }
        a = (ll)(a * a) % MOD; b >>= 1;
    }
    return c;
}

int main(void) {
    freopen("supergcd.in", "r", stdin);
    freopen("supergcd.out", "w", stdout);
    scanf("%d", &n);
    for (int i = 1; i <= n; i++) {
        scanf("%d", &a[i]); lim = min(lim, a[i]);
    }
    mu[1] = 1;
    for (int i = 2; i <= lim; i++) {
        if (!vis[i]) {
            prime[pcnt++] = i; mu[i] = -1;
        }
        for (int j = 0; j < pcnt && (x = prime[j] * i) <= lim; j++) {
            vis[x] = 1;
            if (i % prime[j]) {
                mu[x] = -mu[i];
            } else {
                mu[x] = 0; break;
            }
        }
    }
    for (int i = 1; i <= lim; i++) po[i] = pow(i, i);
    for (int i = 1; i <= lim; i++)
        for (int j = i; j <= lim; j += i)
            f[j] = (f[j] + po[i] * mu[j / i] + MOD) % MOD;
    for (int i = 1; i <= lim; i++) pre[i] = (pre[i - 1] + f[i]) % MOD; //预处理f(G)函数及其前缀和
    for (int i = 1; i <= lim; i = pos + 1) {
        pos = lim;
        for (int j = 1; j <= n; j++) pos = min(pos, a[j] / (a[j] / i));
        res = (pre[pos] - pre[i - 1] + MOD) % MOD;
        for (int j = 1; j <= n; j++) res = (res * (a[j] / i)) % MOD;
        ans = (ans + res) % MOD;
    } //莫比乌斯反演
    cout << ans << endl;
    fclose(stdin); fclose(stdout);
    return 0;
}