原题链接: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 之和。
分析
转化一下题意,列出式子
∑
a
1
=
1
m
∑
a
2
=
1
m
.
.
.
∑
a
n
=
1
m
(
a
1
a
2
a
3
.
.
.
a
n
)
k
[
g
c
d
(
a
1
,
a
2
.
.
.
a
n
)
=
=
d
]
\sum_{a_1=1}^{m}\sum_{a_2=1}^{m}...\sum_{a_n=1}^{m}(a_1a_2a_3...a_n)^k[gcd(a_1,a_2...a_n)==d]
a1=1∑ma2=1∑m...an=1∑m(a1a2a3...an)k[gcd(a1,a2...an)==d]
提出一个d
∑
a
1
=
1
m
/
d
∑
a
2
=
1
m
/
d
.
.
.
∑
a
n
=
1
m
/
d
d
n
k
(
a
1
a
2
a
3
.
.
.
a
n
)
k
[
g
c
d
(
a
1
,
a
2
.
.
.
a
n
)
=
=
1
]
\sum_{a_1=1}^{m/d}\sum_{a_2=1}^{m/d}...\sum_{a_n=1}^{m/d}d^{nk}(a_1a_2a_3...a_n)^k[gcd(a_1,a_2...a_n)==1]
a1=1∑m/da2=1∑m/d...an=1∑m/ddnk(a1a2a3...an)k[gcd(a1,a2...an)==1]
对gcd反演
∑
a
1
=
1
m
/
d
∑
a
2
=
1
m
/
d
.
.
.
∑
a
n
=
1
m
/
d
d
n
k
(
a
1
a
2
a
3
.
.
.
a
n
)
k
∑
a
1
∣
t
,
a
2
∣
t
.
.
μ
(
t
)
\sum_{a_1=1}^{m/d}\sum_{a_2=1}^{m/d}...\sum_{a_n=1}^{m/d}d^{nk}(a_1a_2a_3...a_n)^k\sum_{a_1|t ,a2|t..}\mu(t)
a1=1∑m/da2=1∑m/d...an=1∑m/ddnk(a1a2a3...an)ka1∣t,a2∣t..∑μ(t)
交换枚举顺序
∑
t
=
1
m
/
d
t
n
k
μ
(
t
)
∑
a
1
=
1
m
/
t
d
∑
a
2
=
1
m
/
t
d
.
.
.
∑
a
n
=
1
m
/
t
d
d
n
k
(
a
1
a
2
a
3
.
.
.
a
n
)
k
\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_1a_2a_3...a_n)^k
t=1∑m/dtnkμ(t)a1=1∑m/tda2=1∑m/td...an=1∑m/tddnk(a1a2a3...an)k
观察后面的式子,利用多项式知识可以化简为
∑
t
=
1
m
/
d
t
n
k
μ
(
t
)
d
n
k
(
1
k
+
2
k
+
.
.
.
+
(
m
/
t
d
)
k
)
n
\sum_{t=1}^{m/d}t^{nk}\mu(t)d^{nk}(1^k+2^k+...+(m/td)^k)^n
t=1∑m/dtnkμ(t)dnk(1k+2k+...+(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();
}