【ybt金牌导航8-4-4】【luogu P2155】互质个数 / 沙拉公主的困惑(gcd)(欧拉函数)

互质个数 / 沙拉公主的困惑

题目链接: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=1n![gcd(i,m!)=1](modr)

我们不妨从 n > m n>m n>m 这个方面下手:
那首先 1 ∼ m ! 1\sim m! 1m! 的部分肯定是 φ ( 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! 1m! 的部分循环下去,然后因为是 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!pipi1)=n!pipi1
p i p_i pi m ! m! m! 质因数分解中的每个质数,其实就相当于 1 ∼ m 1\sim m 1m 中的质数)

然后就可以直接搞了(用欧拉筛求质数的时候可以把 m m m 取每个值时的 ∏ p i − 1 p i \prod\dfrac{p_i-1}{p_i} pipi1 算出来)
吗?

其实会有一个问题,就是 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 nr m < r m<r m<r 或者 n ⩾ 2 r n\geqslant 2r n2r(因为 ∏ \prod 下面至多只会出现一次)。

那如果 ⩾ r \geqslant r r 了但是答案不是 0 0 0,那我们计算的时候就直接把 r r r 的位置略过去。
(就是 n ! n! n! 中和 ∏ \prod 中的下面部分,上面还是要留着 r − 1 r-1 r1

然后就可以了。

代码

#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;
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值