「SOL」小A与两位神仙(洛谷)
博客赶工……
# 题面
给定正整数 \(P\)，保证 \(P\) 为某个素数的幂。给定多组正整数 \(x,y\)，保证 \((x,P)=0,(y,P)=0\)，求是否存在正整数 \(a\)，使得 \(x^a\equiv y\pmod P\)。
数据规模：\(P\le10^{18}\)，询问组数不超过 \(2\times10^4\)。
# 解析
先随便什么方法把 \(P\) 给分解成 \(p^k\)，可以得到 \(\varphi(P)=p^{k-1}(p-1)\)，同时可以说明 \(P\) 一定有原根，记为 \(g\)。
于是任何与 \(P\) 互质的数都可以表示成原根的幂，若 \(g^b=a\)，则记 \({\rm ind}\ a=b\)。
根据欧拉定理，可以知道原方程等价于
不妨看看 \(\rm{ind}\ x\) 有什么性质。对所有 \(x\) 在模 \(P\) 意义下的幂定义集合
可以知道 \(x^i=g^{i\cdot{\rm ind}\ x}\)，即 \({\rm ind}\ {x^i}=i\cdot{\rm ind}\ x\)，并且 \(i\) 可以取任意自然数。又因为 \(x^{\varphi(P)}\equiv 1\)，于是可以得到
而 \(|\mathbb{S}|\) 就是 \(x\) 模 \(P\) 意义下的阶 \({\rm ord}\ x\)，于是阶一定是 \(\varphi(P)\) 的因数。
因为 \(y\equiv x^a\)，所以 \({\rm ind}\ y=a\cdot{\rm ind}\ x\)，那么：
于是「存在 \(a\) 使得 \(x^a\equiv y\)」等价于「\({\rm ord}\ y\mid{\rm ord}\ x\)」。
最后一个问题就是，怎么求阶？根据上面的推导以及阶的定义，阶一定是 \(\varphi(P)\) 的因数，而且 \(\forall {\rm ord}\ x\mid a,x^a\equiv 1\)。
所以可以用试除法，设 \({\rm ord}\ x\) 的初值 \(x_0\) 为 \(\varphi(P)\)，将其试除一个 \(\varphi(P)\) 的质因子 \(p_0\)，若 \(p_0\mid 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}\equiv1\pmod p$，则 $p$ 可能是素数；若不成立，则 $p$ 一定不是素数。
二次探测
若 $x^2\equiv 1\pmod p$ 的解有且仅有 $x\equiv\pm1$，则 $p$ 可能是素数，否则一定不是素数。
可以看出基于这两个原理，Miller_rabin 只是一个随机算法，但是它的正确性在 OI 的数据规模内已经够用了~
如何实现 Miller_rabin？首先试除几个（前10个）小质数，排除掉许多合数，加快算法效率。
因为已经试除过 \(2\)，此时的 \(p\) 一定是奇数，可以拆解为 \(p=2^mq\) 的形式，然后进行二次探测。
随机一个数 \(a\)，或者 \(a\) 直接取小质数，通过如下过程实现二次探测：
- 首先检验 \(a^q\) 是否为 \(\pm1\)，若是，则 \(a^q\) 本身就是 \(x^2\equiv1\) 的解，满足二次探测；
- 然后依次检验 \(a^{2q},a^{4q},\dots,a^{2^mq}\)；
- 若检验到 \(a^{2^iq}\) 时，\(a^{2^iq}\) 第一次等于 \(1\)，此时 \(a^{2^{i-1}q}\neq\pm1\)，但是 \(a^{2^{i-1}q}\) 是 \(x^2\equiv1\) 的解，不满足二次探测；
- 若检验到 \(2^{2^iq}\) 时，\(a^{2^iq}\) 第一次等于 \(-1\)，那么 \(a^{2^iq}\) 是 \(x^2\equiv 1\) 的解，满足二次探测；
- 若检验结束后还没有找到 \(x^2\equiv 1\) 的解，即 \(a^{2^iq}\neq\pm1\)，不满足费马小定理的逆定理。
具体可以参照「# 源代码」中的 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\)，\(b\in[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\n",ord(x)%ord(y)? "No":"Yes");
}
return 0;
}
