Description
-
给出非负整数 d d d 和正整数 n n n 的质因数分解式 n = ∏ i = 1 ω p i α i n = \prod_{i = 1}^{\omega} p_i^{\alpha_i} n=∏i=1ωpiαi,请求出
( f d ( n ) = ∑ i = 1 n [ gcd ( i , n ) = 1 ] i d ) m o d ( 1 0 9 + 7 ) \left(f_d(n) = \sum_{i = 1}^n [\gcd(i, n) = 1]\, i^d \right) \bmod (10^9 + 7) (fd(n)=i=1∑n[gcd(i,n)=1]id)mod(109+7) -
0 ≤ d ≤ 100 , 1 ≤ ω ≤ 1 0 3 , 2 ≤ p i ≤ 1 0 9 , 1 ≤ α i ≤ 1 0 9 0\le d\le 100, 1\le \omega\le 10^3, 2\le p_i\le 10^9, 1\le \alpha_i \le 10^9 0≤d≤100,1≤ω≤103,2≤pi≤109,1≤αi≤109。
Solution
f d ( n ) = ∑ i = 1 n [ gcd ( i , n ) = 1 ] i d = ∑ i = 1 n i d ∑ k ∣ gcd ( i , n ) μ ( k ) = ∑ k ∣ n μ ( k ) ∑ i = 1 n [ k ∣ i ] i d = ∑ k ∣ n μ ( k ) k d ∑ i = 1 n k i d \begin{aligned} f_d(n) & = \sum_{i = 1}^n [\gcd(i, n) = 1] \, i^d \\ & = \sum_{i = 1}^n i^d \sum_{k\mid \gcd(i, n)} \mu(k) \\ & = \sum_{k \mid n} \mu(k) \sum_{i = 1}^n [k\mid i] \, i^d \\ & = \sum_{k \mid n} \mu(k) k^d \sum_{i = 1}^{\frac{n}{k}} i^d \end{aligned} fd(n)=i=1∑n[gcd(i,n)=1]id=i=1∑nidk∣gcd(i,n)∑μ(k)=k∣n∑μ(k)i=1∑n[k∣i]id=k∣n∑μ(k)kdi=1∑knid
后面的等幂和直接用伯努利数做就可以了。
我们知道
S
m
(
n
)
=
∑
k
=
0
n
−
1
k
m
=
1
m
+
1
∑
k
=
0
m
(
m
+
1
k
)
B
k
n
m
+
1
−
k
\begin{aligned} S_m(n) & = \sum_{k = 0}^{n - 1} k^m \\ & = \dfrac{1}{m + 1} \sum_{k = 0}^m \dbinom{m + 1}{k} B_k n^{m + 1 - k} \end{aligned}
Sm(n)=k=0∑n−1km=m+11k=0∑m(km+1)Bknm+1−k
设出系数
f
i
=
1
d
+
1
(
d
+
1
i
)
B
i
f_i = \dfrac{1}{d + 1} \dbinom{d + 1}{i} B_i
fi=d+11(id+1)Bi
系数可以
O
(
d
2
)
\Omicron(d^2)
O(d2) 预处理。
则
f
d
(
n
)
=
∑
k
∣
n
μ
(
k
)
k
d
S
d
(
n
k
+
1
)
=
∑
k
∣
n
μ
(
k
)
k
d
∑
i
=
0
d
f
i
(
n
k
+
1
)
d
+
1
−
i
\begin{aligned} f_d(n) & = \sum_{k\mid n} \mu(k) k^d S_d\left(\dfrac{n}{k} + 1 \right) \\ & = \sum_{k\mid n} \mu(k) k^d \sum_{i = 0}^d f_i \left(\dfrac{n}{k} + 1 \right)^{d + 1 - i} \end{aligned}
fd(n)=k∣n∑μ(k)kdSd(kn+1)=k∣n∑μ(k)kdi=0∑dfi(kn+1)d+1−i
你发现这个
(
n
k
+
1
)
m
+
1
−
i
\left(\dfrac{n}{k} + 1 \right)^{m + 1 - i}
(kn+1)m+1−i 没法处理,但如果是
(
n
k
)
m
+
1
−
i
\left(\dfrac{n}{k}\right)^{m + 1 - i}
(kn)m+1−i 就很好了——它可以和前面的
k
d
k^d
kd 相乘。
退回到上一步
f
d
(
n
)
=
∑
k
∣
n
μ
(
k
)
k
d
∑
i
=
1
n
k
i
d
=
∑
k
∣
n
μ
(
k
)
k
d
∑
i
=
1
n
k
−
1
i
d
+
∑
k
∣
n
μ
(
k
)
k
d
(
n
k
)
d
=
∑
k
∣
n
μ
(
k
)
k
d
∑
i
=
1
n
k
−
1
i
d
+
n
d
∑
k
∣
n
μ
(
k
)
=
∑
k
∣
n
μ
(
k
)
k
d
∑
i
=
1
n
k
−
1
i
d
+
n
d
[
n
=
1
]
\begin{aligned} f_d(n) & = \sum_{k\mid n} \mu(k) k^d \sum_{i = 1}^{\frac{n}{k}} i^d \\ & = \sum_{k\mid n} \mu(k) k^d \sum_{i = 1}^{\frac{n}{k} - 1} i^d + \sum_{k\mid n} \mu(k) k^d \left(\dfrac{n}{k}\right)^d \\ & = \sum_{k\mid n} \mu(k) k^d \sum_{i = 1}^{\frac{n}{k} - 1} i^d + n^d \sum_{k\mid n} \mu(k) \\ & = \sum_{k\mid n} \mu(k) k^d \sum_{i = 1}^{\frac{n}{k} - 1} i^d + n^d [n = 1] \\ \end{aligned}
fd(n)=k∣n∑μ(k)kdi=1∑knid=k∣n∑μ(k)kdi=1∑kn−1id+k∣n∑μ(k)kd(kn)d=k∣n∑μ(k)kdi=1∑kn−1id+ndk∣n∑μ(k)=k∣n∑μ(k)kdi=1∑kn−1id+nd[n=1]
观察数据范围
1 ≤ ω ≤ 1 0 3 , 2 ≤ p i ≤ 1 0 9 , 1 ≤ α i ≤ 1 0 9 1\le \omega\le 10^3, 2\le p_i\le 10^9, 1\le \alpha_i \le 10^9 1≤ω≤103,2≤pi≤109,1≤αi≤109。
说明 n ≠ 1 n\ne 1 n=1。
则
f
d
(
n
)
=
∑
k
∣
n
μ
(
k
)
k
d
∑
i
=
1
n
k
−
1
i
d
=
∑
k
∣
n
μ
(
k
)
k
d
S
d
(
n
k
)
=
∑
k
∣
n
μ
(
k
)
k
d
∑
i
=
0
d
f
i
(
n
k
)
d
+
1
−
i
=
∑
k
∣
n
μ
(
k
)
∑
i
=
0
d
f
i
n
d
+
1
−
i
k
i
−
1
=
∑
i
=
0
d
f
i
n
d
+
1
−
i
∑
k
∣
n
μ
(
k
)
k
i
−
1
\begin{aligned} f_d(n) & = \sum_{k\mid n} \mu(k) k^d \sum_{i = 1}^{\frac{n}{k} - 1} i^d \\ & = \sum_{k \mid n} \mu(k) k^d S_d\left(\dfrac{n}{k} \right) \\ & = \sum_{k\mid n} \mu(k) k^d \sum_{i = 0}^d f_i \left(\dfrac{n}{k}\right)^{d + 1 - i} \\ & = \sum_{k\mid n} \mu(k) \sum_{i = 0}^d f_i n^{d + 1 - i} k^{i - 1} \\ & = \sum_{i = 0}^d f_i n^{d + 1 - i} \sum_{k\mid n} \mu(k) k^{i - 1} \end{aligned}
fd(n)=k∣n∑μ(k)kdi=1∑kn−1id=k∣n∑μ(k)kdSd(kn)=k∣n∑μ(k)kdi=0∑dfi(kn)d+1−i=k∣n∑μ(k)i=0∑dfind+1−iki−1=i=0∑dfind+1−ik∣n∑μ(k)ki−1
你会发现
∑
k
∣
n
μ
(
k
)
k
i
−
1
=
[
(
μ
⋅
Id
i
−
1
)
∗
1
]
(
n
)
\sum_{k\mid n} \mu(k) k^{i - 1} = [(\mu\cdot \operatorname{Id}_{i - 1}) * \mathbf{1}](n)
k∣n∑μ(k)ki−1=[(μ⋅Idi−1)∗1](n)
所以这是一个积性函数,但是也不可能线性筛啊。
我们设
∑
k
∣
n
μ
(
k
)
k
i
−
1
=
g
i
(
n
)
\sum_{k\mid n} \mu(k) k^{i - 1} = g_i(n)
k∣n∑μ(k)ki−1=gi(n)
因为
g
i
(
n
)
g_i(n)
gi(n) 是积性函数,所以可以按质因数拆开:
g
i
(
n
)
=
∏
k
=
1
ω
g
i
(
p
k
α
k
)
g_i(n) = \prod_{k = 1}^{\omega} g_i(p_k^{\alpha_k})
gi(n)=k=1∏ωgi(pkαk)
则只用考虑质因数幂的情况。
你会发现
g
i
(
p
α
)
=
∑
k
∣
p
α
μ
(
k
)
k
i
−
1
g_i(p^{\alpha}) = \sum_{k\mid p^{\alpha}} \mu(k) k^{i - 1}
gi(pα)=k∣pα∑μ(k)ki−1
其中只有
k
=
1
k = 1
k=1 和
k
=
p
k = p
k=p 时
μ
(
k
)
\mu(k)
μ(k) 才有贡献。
而
{
μ
(
1
)
1
i
−
1
=
1
μ
(
p
)
p
i
−
1
=
−
p
i
−
1
\begin{cases} \mu(1) 1^{i - 1} = 1 \\ \mu(p) p^{i - 1} = - p^{i - 1} \end{cases}
{μ(1)1i−1=1μ(p)pi−1=−pi−1
所以
g
i
(
p
α
)
=
1
−
p
i
−
1
g_i(p^{\alpha}) = 1 - p^{i - 1}
gi(pα)=1−pi−1
再贴一遍式子
f
d
(
n
)
=
∑
i
=
0
d
f
i
n
d
+
1
−
i
∑
k
∣
n
μ
(
k
)
k
i
−
1
f_d(n) = \sum_{i = 0}^d f_i n^{d + 1 - i} \sum_{k\mid n} \mu(k) k^{i - 1}
fd(n)=i=0∑dfind+1−ik∣n∑μ(k)ki−1
计算
g
i
(
n
)
=
∏
k
=
1
ω
1
−
p
k
i
−
1
g_i(n) = \prod_{k = 1}^{\omega} 1 - p_k^{i - 1}
gi(n)=k=1∏ω1−pki−1
是
O
(
ω
log
i
)
\Omicron(\omega \log i)
O(ωlogi) 的,求个和就是
O
(
ω
d
log
d
)
\Omicron(\omega d \log d)
O(ωdlogd),但其实可以把
p
i
−
1
p^{i - 1}
pi−1 往上递推做到
O
(
d
ω
)
\Omicron(d\omega)
O(dω)。
预处理伯努利数 f i f_i fi 是 O ( d 2 ) \Omicron(d^2) O(d2) 的,而 n d + 1 − i n^{d + 1 - i} nd+1−i 是 O ( d log d ) \Omicron(d\log d) O(dlogd) 的,远不及前面的 O ( d ω ) \Omicron(d\omega) O(dω)。
综上,时间复杂度为 O ( d 2 + d ω ) \Omicron(d^2 + d\omega) O(d2+dω)。
Code
注意预处理要处理到 d + 1 d + 1 d+1。
//18 = 9 + 9 = 18.
#include <iostream>
#include <cstdio>
#define Debug(x) cout << #x << "=" << x << endl
typedef long long ll;
using namespace std;
const int MOD = 1e9 + 7;
const int MAXD = 105;
int add(int a, int b) {return (a + b) % MOD;}
int sub(int a, int b) {return (a - b + MOD) % MOD;}
int mul(int a, int b) {return (ll)a * b % MOD;}
int qpow(int a, int b) {int base = a, ans = 1; while (b) {if (b & 1) ans = mul(ans, base); base = mul(base, base); b >>= 1;} return ans;}
int GetInv(int a) {return qpow(a, MOD - 2);}
struct Bernoulli
{
int fac[MAXD], fac_inv[MAXD], inv[MAXD];
void pre(int d)
{
fac[0] = fac_inv[0] = inv[1] = 1;
for (int i = 1; i <= d; i++)
{
fac[i] = mul(fac[i - 1], i);
}
fac_inv[d] = GetInv(fac[d]);
for (int i = d - 1; i >= 1; i--)
{
fac_inv[i] = mul(fac_inv[i + 1], i + 1);
inv[i + 1] = mul(fac_inv[i + 1], fac[i]);
}
}
int C(int n, int m)
{
return mul(mul(fac[n], fac_inv[n - m]), fac_inv[m]);
}
int B[MAXD], f[MAXD];
void cal(int d)
{
B[0] = 1;
f[0] = inv[d + 1];
for (int m = 1; m <= d; m++)
{
int sum = 0;
for (int k = 0; k < m; k++)
{
sum = add(sum, mul(C(m + 1, k), B[k]));
}
B[m] = mul(MOD - sum, inv[m + 1]);
f[m] = mul(mul(inv[d + 1], C(d + 1, m)), B[m]);
}
}
}Ber;
const int MAXW = 1005;
int p[MAXW], a[MAXW], xsy[MAXW];
int main()
{
int d, w;
scanf("%d%d", &d, &w);
Ber.pre(d + 1), Ber.cal(d);
int n = 1;
for (int i = 1; i <= w; i++)
{
scanf("%d%d", p + i, a + i);
n = mul(n, qpow(p[i], a[i]));
xsy[i] = GetInv(p[i]);
}
int ans = 0;
for (int i = 0; i <= d; i++)
{
int pro = 1;
for (int j = 1; j <= w; j++)
{
pro = mul(pro, sub(1, xsy[j]));
xsy[j] = mul(xsy[j], p[j]);
}
ans = add(ans, mul(mul(Ber.f[i], qpow(n, d + 1 - i)), pro));
}
printf("%d\n", ans);
return 0;
}