TJOI 2017 (Pollard_rho,MR素数测试）

WedsonLin

于 2020-12-13 17:18:00 发布

阅读量201

点赞数

CC 4.0 BY-SA版权

分类专栏：数论 ACM 文章标签： acm acm竞赛

本文链接：https://blog.youkuaiyun.com/WedsonLin/article/details/111129556

ACM 同时被 2 个专栏收录

41 篇文章

订阅专栏

数论

25 篇文章

订阅专栏

博客详细介绍了如何在TJOI 2017竞赛中使用Pollard_rho算法和Miller-Rabin素数测试解决求分数AB模M逆元的问题。通过分析质因子分解，简化求逆过程，避免TLE。最终代码使用了__int128数据类型来处理大数。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

TJOI 2017 (Pollard_rho,MR素数测试）

题目链接

题意是求分数 $\frac{B}{A}$ 在模M下的逆元，没有则输出-1。

本题的两大核心是：Pollard_rho算法，和MR素数测试算法。

$\frac{B}{A}$ 在模M下有逆元当且仅当 $g c d (A, M) = 1$ ，因为分数不是最简的，所以只要分数可以化简到 $\frac{B'}{A'}$ ，其中 $g c d (A^{'}, M) = 1$ 即有逆元。

我一开始的思路是先将 $\frac{B}{A}$ 化到最简分数 $\frac{B'}{A'}$ ，再判断A’是否有逆元。由于B’和A‘超出了2⁶⁴的范围，所以不能直接求 $g c d$ 以及约分。因为 $A=a_1\cdot a_2\cdots a_m，B=b_1\cdot b_2\cdots b_m$ ，所以可以通过依次求a_i,b_i的质因子分解式，将A,B的质因子分解式求出，在通过A,B的质因子分解式去化简 $\frac{B}{A}$ ，再求M的质因子分解式，这样就可以判断是否 $g c d (A^{'}, M) = 1$ 了。

思路是简单明了的，答案能也保证正确，但是TLE。

在参考了别人的代码后，我发现并不用求A,B的质因子分解式，甚至也不用求M的质因子分解式，只需要求M的质因子集合即可。其实并不需要将 $\frac{B}{A}$ 化到最简，只需要将 $\frac{B}{A}$ 化简到 $\frac{B'}{A'}$ 使得 $g c d (A^{'}, M) = 1$ 即可。

在得到M的质因子集合 $s=\{p_1,p_2,\cdots p_n\}$ 之后，想要让 $g c d (A^{'}, M) = 1$ ，那么A’必须不能包含集合s里面的质因子。因为 $A=a_1\cdot a_2\cdots a_m$ ，所以我们依次检查a_i是否能被p_j整除，若能整除，则一定有 $a_i=a^{'}_i*p_j^{e_j}$ ,保留 $a^{'}$ 并统计 $e^j$ 。得到的 $A'=a^{'}_1\cdot a^{'}_2\cdots a^{'}_n$ 满足于M互质。
因为我们需要保证变化后的分数是等价的，所以每当A除以一个p_j时，我们还需判断B是否也能除以一个p_j。若不能，则说明 $\frac{B}{A}$ 化简后的 $\frac{B'}{A'}$ 中的A’中必含有集合s里面至少一个素因子。所以我们还需统计 $B=b_1\cdot b_2\cdots b_m，b_i=b^{'}_i*p_j^{e_j}$ 中的e_j。

可以用一个数组e，初始化e[i]=0：
对于 $a_i=a^{'}_i*p_j^{e_j}$ ,e[j]-=e_j。
对于 $b_i=b^{'}_i*p_j^{e_j}$ ，e[j]+=e_j。

若存在e[j]<0，则说明 $\frac{B}{A}$ 并不存在逆元。
若对于所有e[j]>=0，则一定存在逆元。

注意，最后得到的 $A'=a^{'}_1\cdot a^{'}_2\cdots a^{'}_n$ , $B'=b^{'}_1\cdot b^{'}_2\cdots b^{'}_n$ 组成的 $\frac{B'}{A'}$ 并不等价于原来的分数，因为可能存在e[i]>0,说明了B在约去了A中包含的质因子p_i之外还多约了 $p_i^{{e[i]}}$ 。所以需将多约去的返还给B。

（不明白题目保证1<M,ai,bi<2*10^8，但是64bits的long long 过不了？）
所以本题实现使用了__int128。

代码如下：

#include<iostream>
#include<cstdio>
#include<utility>
#include<unordered_map>
#include<vector>
#include<algorithm>
#include<cstdlib>
#include<ctime>
using namespace std;
typedef long long ll;
const int max_n=1e4+5; 
ll b[max_n];
ll a[25][max_n];
ll c[55];
ll M[55];
ll aa[max_n],bb[max_n];
ll n,m,k;
vector<ll> f;
ll gcd(ll a,ll b)
{
	ll t;
	while(b)
	{
		t=a;
		a=b;
		b=t%b;
	}
	return a;
}
ll qpow(ll a,ll n,ll mod)
{
	ll ans=1;
	while(n)
	{
		if(n&1)ans=(__int128)ans*a%mod;
		a=(__int128)a*a%mod;
		n>>=1;
	}
	return ans;
}
bool mr_test(ll n)
{
	if(n<2)return false;
	else if(n==2)return true;
	else if(n%2==0)return false;
	ll a=n-1,b=0;
	while(a%2==0)a/=2,b++;
	for(int i=0;i<10;i++)
	{
		ll x=rand()%(n-2)+2;
		ll y=qpow(x,a,n);
		for(int j=0;j<b;j++)
		{
			x=(__int128)y*y%n;
			if(x==1&&y!=1&&y!=n-1)return false;
			y=x;
		}
		if(y!=1)return false;
	}
	return true;
}
void pollard_rho(ll n)
{
	if(n==1)return ;
	if(mr_test(n)){
		f.push_back(n);
		return ;
	}
	ll c=rand()%(n-1)+1;
	ll s=0,t=0;
	ll v=1;
	ll goal,step=1;
	for(goal=1;;goal<<=1,s=t,v=1)
	{
		for(;step<=goal;step++)
		{
			t=((__int128)t*t%n+c)%n;
			v=(__int128)v*abs(t-s)%n;
			if(step%127==0)
			{
				ll d=gcd(v,n);
				if(d>1){
					while(n%d==0)n/=d;
					pollard_rho(n);
					pollard_rho(d);
					return ;
				}
			}
		}
		ll d=gcd(v,n);
			if(d>1){
				while(n%d==0)n/=d;
				pollard_rho(n);
				pollard_rho(d);
				return ;
			}		
	}
}
void solve(void)
{
	for(int i=1;i<=k;i++,f.clear())
	{
		pollard_rho(M[i]);
		sort(f.begin(),f.end());
		int cnt=unique(f.begin(),f.end())-f.begin();
		vector<ll> e(cnt);
		for(int j=1;j<=m;j++)aa[j]=a[c[i]][j],bb[j]=b[j];
		bool flag=false;
		for(int j=0;j<cnt;j++)
		{
			for(int l=1;l<=m;l++)while(bb[l]%f[j]==0)bb[l]/=f[j],e[j]++;
			for(int l=1;l<=m;l++)while(aa[l]%f[j]==0)aa[l]/=f[j],e[j]--;
			if(e[j]<0){
				flag=true;
				break;
			}
		}
		if(flag){
			printf("-1\n");
			continue;
		}
		ll ansa=1,ansb=1,phi=M[i];
		for(int j=0;j<cnt;j++)phi-=phi/f[j];
		for(int j=1;j<=m;j++)ansa=(__int128)ansa*aa[j]%M[i];
		for(int j=1;j<=m;j++)ansb=(__int128)ansb*bb[j]%M[i];
		for(int j=0;j<cnt;j++)ansb=(__int128)ansb*qpow(f[j],e[j],M[i])%M[i];//把除多的补回来 
		ll rev=qpow(ansa,phi-1,M[i]);
		ll ans=(__int128)ansb*rev%M[i];
		printf("%lld\n",ans);
	}
	return ;
}
int main(void)
{
	srand(time(NULL));
	scanf("%lld%lld%lld",&n,&m,&k);
	for(int i=1;i<=m;i++)
	scanf("%lld",&b[i]);
	for(int i=1;i<=n;i++)
	for(int j=1;j<=m;j++)
	scanf("%lld",&a[i][j]);
	for(int i=1;i<=k;i++)
	scanf("%lld %lld",&c[i],&M[i]);
	solve();
}