2019 ICPC Asia Yinchuan Regional D - Easy Problem 莫比乌斯反演 + 欧拉降幂

原题链接：https://nanti.jisuanke.com/t/42384

题意

给定一个序列，我们称一个序列最大公因数为d，元素个数为n，大小不超过m的序列为（n, m, k）三元组，设序列a1,a2,a3…an为q，求f（q, k）表示所有满足条件的（a1, a2, a3…an）^k 之和。

分析

转化一下题意，列出式子
$$\sum _{a_1=1}^m\sum _{a_2=1}^m...\sum _{a_n=1}^m(a_1{a}_2{a}_3...a_n{)}^k[gcd(a_1,a_2...a_n)==d]$$
提出一个d
$$\sum _{a_1=1}^{m/d}\sum _{a_2=1}^{m/d}...\sum _{a_n=1}^{m/d}{d}^{nk}(a_1{a}_2{a}_3...a_n{)}^k[gcd(a_1,a_2...a_n)==1]$$
对gcd反演
$$\sum _{a_1=1}^{m/d}\sum _{a_2=1}^{m/d}...\sum _{a_n=1}^{m/d}{d}^{nk}(a_1{a}_2{a}_3...a_n{)}^k\sum _{a_1|t,a2|t..}\mu (t)$$
交换枚举顺序
$$\sum _{t=1}^{m/d}{t}^{nk}\mu (t)\sum _{a_1=1}^{m/td}\sum _{a_2=1}^{m/td}...\sum _{a_n=1}^{m/td}{d}^{nk}(a_1{a}_2{a}_3...a_n{)}^k$$
观察后面的式子，利用多项式知识可以化简为

$$\sum _{t=1}^{m/d}{t}^{nk}\mu (t)d^{nk}(1^k+2^k+...+(m/td{)}^k{)}^n$$
由于n的范围很大，要用到欧拉降幂。
后面的式子可以O（n）预处理一遍，然后最后枚举也是O（n）就可以了。

Code

#include <bits/stdc++.h>
using namespace std;
//#define ACM_LOCAL
#define re register
#define fi first
#define se second
const int N = 1e6 + 10;
const int M = 1e6 + 10;
const int INF = 0x3f3f3f3f;
const double eps = 1e-4;
const int MOD = 59964251;
typedef long long ll;
typedef pair<int, int> PII;
ll prime[N], mu[N];
ll phi[N], ff[N], k, K;
bool is_prime[N];
ll ksm(ll a, ll b) {
    ll res = 1, base = a;
    while (b) {
        if (b & 1) res = res * base % MOD;
        base = base * base % MOD;
        b >>= 1;
    }
    return res;
}
void init(int n) {
    memset(is_prime, true, sizeof is_prime);
    mu[1] = 1; phi[1] = 1;
    for (re int i = 2; i < n; ++i) {
        if (is_prime[i]) prime[++k] = i, mu[i] = -1, phi[i] = i-1;
        for (re int j = 1; j <= k && i * prime[j] < n; ++j) {
            is_prime[i * prime[j]] = false;
            if (i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            } else {
                mu[i * prime[j]] = -mu[i];
                phi[i * prime[j]] = phi[i] * (prime[j] - 1);
            }
        }
    }
}
ll get_phi(ll n) {
    ll m = (int)sqrt(n + 0.5);
    ll ans = n;
    for (ll i = 2; i <= m; ++ i) {
        if (n % i == 0) {
            ans = ans / i * (i - 1);
            while(n % i == 0) n /= i;
        }
    }
    if (n > 1) ans = ans / n * (n - 1);
    return ans;
}
ll mod(ll a, ll mm) {return a >= mm ? a % mm + mm : a;}
ll mm, _n;
ll calc(ll m) {
    ll ans = 0;
    ff[0] = 0;
    for (int i = 1; i <= m; ++i) {
        ff[i] = (ff[i-1] + ksm(i, K)) % MOD;
    }
    for (int i = 1; i <= m; i++) {
        ans = (ans + mu[i] * ksm(i, mod(K * _n, mm)) % MOD * ksm(ff[m/i], _n) % MOD + MOD) % MOD;
    }
    return ans;
}
void solve() {
    mm = get_phi(MOD);
    init(1e5 + 5);
    int T; cin >> T; while (T--) {
        string n; ll m, d;
        cin >> n >> m >> d >> K;
        _n = 0; int flag = 0;
        for (int i = 0; i < n.size(); i++) {
            _n = _n * 10 + n[i] - '0';
            if (_n >= mm) {
                _n %= mm;
                flag = 1;
            }
        }
        if (flag) _n += mm, _n %= mm;
        printf("%lld\n", ksm(d, mod(_n * K, mm)) * calc(m / d) % MOD);
    }
}
signed main() {
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
#ifdef ACM_LOCAL
    freopen("input", "r", stdin);
    freopen("output", "w", stdout);
#endif
    solve();
}