博客赶工……
# 题面
给定正整数 (P)，保证 (P) 为某个素数的幂。给定多组正整数 (x,y)，保证 ((x,P)=0,(y,P)=0)，求是否存在正整数 (a)，使得 (x^aequiv ypmod P)。
数据规模：(Ple10^{18})，询问组数不超过 (2 imes10^4)。
# 解析
先随便什么方法把 (P) 给分解成 (p^k)，可以得到 (varphi(P)=p^{k-1}(p-1))，同时可以说明 (P) 一定有原根，记为 (g)。
于是任何与 (P) 互质的数都可以表示成原根的幂，若 (g^b=a)，则记 ({ m ind} a=b)。
根据欧拉定理，可以知道原方程等价于
不妨看看 ( m{ind} x) 有什么性质。对所有 (x) 在模 (P) 意义下的幂定义集合
可以知道 (x^i=g^{icdot{ m ind} x})，即 ({ m ind} {x^i}=icdot{ m ind} x)，并且 (i) 可以取任意自然数。又因为 (x^{varphi(P)}equiv 1)，于是可以得到
而 (|mathbb{S}|) 就是 (x) 模 (P) 意义下的阶 ({ m ord} x)，于是阶一定是 (varphi(P)) 的因数。
因为 (yequiv x^a)，所以 ({ m ind} y=acdot{ m ind} x)，那么：
于是「存在 (a) 使得 (x^aequiv y)」等价于「({ m ord} ymid{ m ord} x)」。
最后一个问题就是，怎么求阶？根据上面的推导以及阶的定义，阶一定是 (varphi(P)) 的因数，而且 (forall { m ord} xmid a,x^aequiv 1)。
所以可以用试除法，设 ({ m ord} x) 的初值 (x_0) 为 (varphi(P))，将其试除一个 (varphi(P)) 的质因子 (p_0)，若 (p_0mid x_0) 且 (x^{frac{x_0}{p_0}}equiv 1)，则 (x_0:=frac{x_0}{p_0})。
至于怎么找到 (varphi(P)) 的质因子，由于数很大，可以使用 Pollard_rho。
复习笔记
# Miller_rabin 素性测试
检测原理基于「费马定理的逆定理」以及「二次探测」。
费马定理的逆定理
若 $a^{p-1}equiv1pmod p$，则 $p$ 可能是素数；若不成立，则 $p$ 一定不是素数。
二次探测
若 $x^2equiv 1pmod p$ 的解有且仅有 $xequivpm1$，则 $p$ 可能是素数，否则一定不是素数。
可以看出基于这两个原理，Miller_rabin 只是一个随机算法，但是它的正确性在 OI 的数据规模内已经够用了~
如何实现 Miller_rabin？首先试除几个（前10个）小质数，排除掉许多合数，加快算法效率。
因为已经试除过 (2)，此时的 (p) 一定是奇数，可以拆解为 (p=2^mq) 的形式，然后进行二次探测。
随机一个数 (a)，或者 (a) 直接取小质数，通过如下过程实现二次探测：
- 首先检验 (a^q) 是否为 (pm1)，若是，则 (a^q) 本身就是 (x^2equiv1) 的解，满足二次探测；
- 然后依次检验 (a^{2q},a^{4q},dots,a^{2^mq})；
- 若检验到 (a^{2^iq}) 时，(a^{2^iq}) 第一次等于 (1)，此时 (a^{2^{i-1}q} eqpm1)，但是 (a^{2^{i-1}q}) 是 (x^2equiv1) 的解，不满足二次探测；
- 若检验到 (2^{2^iq}) 时，(a^{2^iq}) 第一次等于 (-1)，那么 (a^{2^iq}) 是 (x^2equiv 1) 的解，满足二次探测；
- 若检验结束后还没有找到 (x^2equiv 1) 的解，即 (a^{2^iq} eqpm1)，不满足费马小定理的逆定理。
具体可以参照「# 源代码」中的 miller_rabin 函数。
# Pollard_rho 因数分解
用于找到一个非 (1) 正整数 (P) 的大于 (1) 小于 (P) 的因子，可以辅助分解质因数。
Pollard_rho 也为一个随机算法，进行一次 Pollard_rho 有可能无法找出这样的因子。
Pollard_rho 的概率保证主要有以下两点：
- 随机 (a,b) 计算 (|a-b|)，而不是直接随机一个数 (a)（生日悖论，但是我也不是很懂这个）；
- 并不直接判断 (|a-b|) 是 (P) 的因数，而是判断 (|a-b|) 是否与 (P) 互质，若不互质，则它们的最大公约数就是 (P) 的一个大于 (1) 的因子，这样符合条件的 (|a-b|) 就更多了。
所以实现时，我们采用伪随机数的方式——令伪随机函数 (f(x)=x^2+b)，(bin[1,P))。
先随机两个变量 (x_1=x_2)，然后每次以不同的倍数对 (x_1,x_2) 进行伪随机化，例如 (x_1:=f(x_1),x_2:=f(f(x_2)))。
可以发现这样 (x_1,x_2) 在模 (P) 意义下一定会产生循环，并且在循环内可能找不到 (|x_1-x_2|) 与 (P) 不互质。这就是 Pollard_rho 是随机算法的原因。
产生循环时 (x_1=x_2)，此时 ((|x_1-x_2|,P)=P)，自动退出循环。
另外，当然不能对一个质数进行 Pollard_rho 算法，因此先用 Miller_rabin 把 (P) 是质数的情况判掉。
具体可以参照「# 源代码」中的 pollard_rho 函数。
# 源代码
点击展开/折叠代码
/*Lucky_Glass*/
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
template<class T>inline T rin(T &r){
int b=1,c=getchar();r=0;
while(c<'0' || '9'<c) b=c=='-'?-1:b,c=getchar();
while('0'<=c && c<='9') r=(r<<1)+(r<<3)+(c^'0'),c=getchar();
return r*=b;
}
typedef long long llong;
#define con(type) const type &
namespace PRIME{
llong ina_gcd(con(llong)a,con(llong)b){return b?ina_gcd(b,a%b):a;}
llong ina_abs(con(llong)a){return a<0?-a:a;}
llong mul(con(llong)a,con(llong)b,con(llong)mod){return llong((__int128)a*b%mod);}
llong ina_pow(llong a,llong b,con(llong)mod){
llong r=1;
while(b){
if(b&1) r=mul(r,a,mod);
a=mul(a,a,mod),b>>=1;
}
return r;
}
bool miller_rabin(con(llong)x,con(llong)a,llong b){
llong tmp=ina_pow(a,b,x);
if(tmp==x-1 || tmp==1) return true;
while(b!=x-1){
tmp=mul(tmp,tmp,x),b<<=1;
if(tmp==x-1) return true;
if(tmp==1) return false;
}
return false;
}
const int PRM[]={2,3,5,7,11,13,17,19,61,24251};
bool primeTest(con(llong)x){
for(int i=0;i<10;i++)
if(x%PRM[i]==0)
return x==PRM[i];
llong tmp=x-1;
while(!(tmp&1)) tmp>>=1;
for(int i=0;i<10;i++)
if(!miller_rabin(x,PRM[i],tmp))
return false;
return true;
}
llong pollard_rho(con(llong)x){
for(int i=0;i<10;i++)
if(x%PRM[i]==0)
return PRM[i];
llong x1=rand()%(x-2)+2,x2=x1,tmp=rand()%(x-1)+1,d=1;
while(d==1){
x1=(mul(x1,x1,x)+tmp)%x;
x2=(mul(x2,x2,x)+tmp)%x;
x2=(mul(x2,x2,x)+tmp)%x;
d=ina_gcd(ina_abs(x1-x2),x);
}
return d;
}
void divideNum(con(llong)x,llong *&res){
if(x==1) return;
if(primeTest(x)){
*res=x,res++;
return;
}
llong dv=pollard_rho(x);
divideNum(dv,res),divideNum(x/dv,res);
}
}
llong m,phim,dv_phim[100],ena_m[100];
int ncas,ndv_phim,nena_m;
llong eta_pow(con(llong)x,con(int)k){
llong ret=1;
for(int i=1;i<=k;i++){
__int128 tmp=(__int128)x*ret;
if(tmp>m) return m+1;
ret=(llong)tmp;
}
return ret;
}
llong kthRoot(con(int)k){
llong tmp=pow(m,1.0/k)+2;
while(eta_pow(tmp,k)>m) tmp--;
return tmp;
}
llong divideM(){
llong tmp;
for(int i=60;i;i--)
if(tmp=kthRoot(i)){
if(eta_pow(tmp,i)==m)
return tmp;
}
return m;
}
void init(){
llong pm=divideM();
phim=pm-1;
for(llong it=pm;it<m;it*=pm)
phim*=pm;
llong *tmp=dv_phim;
PRIME::divideNum(phim,tmp);
ndv_phim=int(tmp-dv_phim);
sort(dv_phim,tmp);
for(int i=0;i<ndv_phim;i++){
int j=i;
while(j+1<ndv_phim && dv_phim[j+1]==dv_phim[i]) j++;
ena_m[++nena_m]=dv_phim[i];
i=j;
}
}
llong ord(con(llong)x){
llong ret=phim;
for(int i=1;i<=nena_m;i++)
while(ret%ena_m[i]==0 && PRIME::ina_pow(x,ret/ena_m[i],m)==1)
ret/=ena_m[i];
return ret;
}
int main(){
// freopen("input.in","r",stdin);
rin(m),rin(ncas);
init();
while(ncas--){
llong x,y;rin(x),rin(y);
printf("%s
",ord(x)%ord(y)? "No":"Yes");
}
return 0;
}