博客赶工……
# 题面
给定正整数 $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 函数。
# 源代码
点击展开/折叠代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
|
#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(){
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;
}
|
THE END
Thanks for reading!
去听从内心的旨意 去步步为营
绝望中爆发无限潜力 为弱肉强食
——《陷落之序》By 著小生zoki
> Link 陷落之序-Bilibili