互质个数 / 沙拉公主的困惑
题目链接:ybt金牌导航8-4-4 / luogu P2155
题目大意
要你求 sum i=1~n![gcd(i,m!)=1]%r。
其中 r 是质数,n>=m。
很多组询问。
思路
重新写式子:
∑
i
=
1
n
!
[
g
c
d
(
i
,
m
!
)
=
1
]
(
m
o
d
r
)
\sum\limits_{i=1}^{n!}[gcd(i,m!)=1]\pmod r
i=1∑n![gcd(i,m!)=1](modr)
我们不妨从
n
>
m
n>m
n>m 这个方面下手:
那首先
1
∼
m
!
1\sim m!
1∼m! 的部分肯定是
φ
(
m
!
)
\varphi(m!)
φ(m!)。
那超出的部分呢?
我们再考虑一下会发现
m
!
∣
n
!
m!|n!
m!∣n!。
然后再根据
gcd
(
a
+
b
,
b
)
=
gcd
(
a
,
b
)
\gcd(a+b,b)=\gcd(a,b)
gcd(a+b,b)=gcd(a,b) 这个东西。
你就会发现超出的部分就是不停的像
1
∼
m
!
1\sim m!
1∼m! 的部分循环下去,然后因为是
m
!
∣
n
!
m!|n!
m!∣n!,所以最后一次一定是循环完了的!
所以就是循环了
n
!
m
!
\dfrac{n!}{m!}
m!n! 次,那答案就是
n
!
m
!
φ
(
m
!
)
\dfrac{n!}{m!}\varphi(m!)
m!n!φ(m!)。
然后直接求好像不太行,因为它是很多组询问,所以我们考虑把
φ
\varphi
φ 用欧拉筛求的原理分解出来。
然后就是:
n
!
m
!
φ
(
m
!
)
=
n
!
m
!
(
m
!
∏
p
i
−
1
p
i
)
=
n
!
∏
p
i
−
1
p
i
\dfrac{n!}{m!}\varphi(m!)=\dfrac{n!}{m!}(m!\prod\dfrac{p_i-1}{p_i})=n!\prod\dfrac{p_i-1}{p_i}
m!n!φ(m!)=m!n!(m!∏pipi−1)=n!∏pipi−1
(
p
i
p_i
pi 是
m
!
m!
m! 质因数分解中的每个质数,其实就相当于
1
∼
m
1\sim m
1∼m 中的质数)
然后就可以直接搞了(用欧拉筛求质数的时候可以把
m
m
m 取每个值时的
∏
p
i
−
1
p
i
\prod\dfrac{p_i-1}{p_i}
∏pipi−1 算出来)
吗?
其实会有一个问题,就是
n
,
m
n,m
n,m 时可能
⩾
r
\geqslant r
⩾r 的。
那这个时候不一定是
0
0
0,因为可能
n
!
n!
n! 会与
∏
\prod
∏ 中的下面的
p
i
p_i
pi 消掉。
所以准确来说无解的情况时
n
⩾
r
n\geqslant r
n⩾r 时
m
<
r
m<r
m<r 或者
n
⩾
2
r
n\geqslant 2r
n⩾2r(因为
∏
\prod
∏ 下面至多只会出现一次)。
那如果
⩾
r
\geqslant r
⩾r 了但是答案不是
0
0
0,那我们计算的时候就直接把
r
r
r 的位置略过去。
(就是
n
!
n!
n! 中和
∏
\prod
∏ 中的下面部分,上面还是要留着
r
−
1
r-1
r−1)
然后就可以了。
代码
#include<cstdio>
#define ll long long
using namespace std;
ll t, mo, n, m, jc[10000001], inv[10000001];
ll prime[1000001];
bool np[10000001];
ll ksm(ll x, ll y) {
ll re = 1;
while (y) {
if (y & 1) re = re * x % mo;
x = x * x % mo;
y >>= 1;
}
return re;
}
int main() {
scanf("%lld %lld", &t, &mo);
jc[0] = 1;
for (int i = 1; i <= 10000000; i++)
if (i != mo) jc[i] = jc[i - 1] * i % mo;
else jc[i] = jc[i - 1];//忽略掉 mo
inv[0] = inv[1] = 1;
for (int i = 2; i <= 10000000; i++) {
if (!np[i]) {
if (i != mo) inv[i] = inv[i - 1] * ksm(i, mo - 2) % mo * (i - 1) % mo;
else inv[i] = inv[i - 1] * (i - 1) % mo;//下面忽略掉,上面还要
prime[++prime[0]] = i;
}
else inv[i] = inv[i - 1];
for (int j = 1; j <= prime[0] && i * prime[j] <= 10000000; j++) {
np[i * prime[j]] = 1;
if (i % prime[j] == 0) break;
}
}
while (t--) {
scanf("%d %d", &n, &m);
if (n >= mo) {
if (m < mo || n >= mo * 2) {//直接判无解情况
printf("0\n");
continue;
}
}
printf("%lld\n", jc[n] * inv[m] % mo);
}
return 0;
}