数论ex
数论ex
数学学得太差了补补知识点or复习
Miller-Rabin 和 Pollard Rho
Miller-Rabin
SRE实战 互联网时代守护先锋，助力企业售后服务体系运筹帷幄！一键直达领取阿里云限量特价优惠。前置知识：
费马小定理
\[ a^{p-1}\equiv 1\pmod p,p \ is \ prime \]二次探测（mod奇素数下1的二次剩余）
\[ x^2\equiv 1\pmod p\Rightarrow x=1 \ or \ p-1 \]
如果不是 \(\bmod\) 奇素数，二次剩余可能是更多的值
如果把费马小定理反过来用来检测一个数是否是素数，虽然是错的，但是反例比较少，如果配合二次探测进行测试并取多个a，可以把1e18内的所有数是否是质数判断出来。
具体的，取\(a\)为前\(12\)个质数\((2,3,5,7,11,13,17,19,23,29,31,37)\)
然后用费马小定理进行测试，如果\(a^{p-1}\equiv 1\pmod p\)，就根据二次探测检验是否有\(a^{(p-1)/2}\equiv 1\pmod p\)，如果值为\(1\)，就递归\((p-1)/4\)处理，如果值为\(p-1\)，无法向下递归直接返回true，否则返回false
Code:
const int pri[]={2,3,5,7,11,13,17,19,23,29,31,37}; bool ck(ll a,ll k,ll p) { ll x=qp(a,k,p);//a^k mod p if(x!=1&&x!=p-1) return false; if(x==p-1) return true; if(k&1) return true; return ck(a,k>>1,p); } bool Miller_Rabin(ll p) { if(p==1) return false; for(int i=0;i<12;i++) if(p%pri[i]==0) return p==pri[i]; for(int i=0;i<12;i++) if(!ck(pri[i],p-1,p)) return false; return true; }这样做的复杂度是\(12\log^2 n\)，其中因子2的个数是一个\(\log\)，里面每次还做了一次快速幂，如果里面用了龟速乘就更完蛋了。
考虑从底向上做，就是说，先求出\(x-1\)除2到不能除的最底层的\(x\)，这样每次向上一层直接就是\(x^2\)了
考虑在某一层\(x=1\)，那么如果之前的一层的\(x'\not=p-1\)的话，就不合法了，或者在最上面一层\(x\)仍然不为\(1\)，也是不合法的
这样就优化掉了一个\(\log\)
Code:
const int pri[]={2,3,5,7,11,13,17,19,23,29,31,37}; bool Miller_Rabin(ll p) { if(p==1) return false; for(int i=0;i<12;i++) if(p%pri[i]==0) return p==pri[i]; ll res=p-1;int k=0; while(!(res&1)) res>>=1,++k; for(int i=0;i<12;i++) { ll x=qp(pri[i],res,p); for(int j=0;j<k&&x>1;j++) { ll y=(ll)x*x%p; if(y==1&&x!=p-1) return false; x=y; } if(x!=1) return false; } return true; }Pollard Rho
我再胡乱说一个自己瞎编的假装是给自己看的，说一下流程把
首先是不加倍增优化的
考虑在\(n\)范围内rand两个数\(x,y\)，如果\(\gcd(x+n-y\bmod n,n)\not=1,n\)，那么\(|x-y|\)就是一个\(n\)的约数。
现在搞一个近似随机函数\(f\)，保证\(f\)在\(\bmod n\)近似均匀随机
令\(f_m(x)=f(f_{m-1}(x))\)，因为\(f\)取值有限，所以一定会出现循环，考虑在出现循环的时候如果我们还没有找到约数，我们就退出。
具体的，初始\(x=y=rand()\)，然后\(x\)一步一步跳，\(x'=f(x)\)，\(y\)两步两步跳\(y'=f(f(y))\)，如果\(x\)与\(y\)又跳到相等了，就说明找到了环，可以证明这个环最多被走一次就找到了。
发现\(f(x)=x^2+c\)的时候效果好，\(c\)是随机正整数。
于是我们每次\(rand()\)一个\(x,y\)和\(c\)去找一个约数，然后分解\(n\)，递归子问题继续找，直到Miller Rabin检测\(n\)为素数。
Code:
ll Find(ll n) { for(int i=0;i<12;i++) if(n%pri[i]==0) return pri[i]; ll x=rand(),y=x,c=rand(); int lim=1e5; do { ll g=gcd((x-y+n)%n,n); if(g!=1&&g!=n) return g; x=f(x,c,n),y=f(f(y,c,n),c,n); }while(x!=y&&lim--); return -1; } void Pollard_Rho(ll n,ll &a) { if(n<a) return; while(!(n&1)) n>>=1;a=2; if(n==1||Miller_Rabin(n)) {a=a>n?a:n;return;} ll d=Find(n); while(d==-1) d=Find(n); if(d<n/d) d=n/d; Pollard_Rho(d,a),Pollard_Rho(n/d,a); }这样的复杂度是\(O(n^{\frac{1}{4}}\log n)\)的
考虑去优化掉\(\log\)
注意到一个事实，\(\gcd(x,n)|\gcd(xy\bmod n,n)\)
于是我们可以倍增\(y\)的步数，然后在每次倍增中，每\(w\)步（\(w\)可取\(127,512\)等值）取一次\(\gcd\)，如果\(\gcd\)值为\(1\)，说明这\(\gcd\)步都没有值，继续向后找，如果\(\gcd\)值为\(n\)，说明可能得到答案，就重新对这\(w\)求一遍，中间如果遇到\(\gcd\not=1,n\)就返回。但这样最后可能返回的值为\(n\)，我们就重新rand()参数继续跑，如果都不是，说明得到了答案，直接返回就可以了。
这样的复杂度是真实\(O(n^{\frac{1}{4}})\)的
Code:
ll Find(ll n) { for(int i=0;i<12;i++) if(n%pri[i]==0) return pri[i]; ll x,y=rand(),c=rand(); int w=1<<9; for(int l=1;;l<<=1) { x=y; for(int i=0;i<l;i++) y=f(y,c,n); for(int i=0;i<l;i+=w) { int le=min(l-i,w); ll g=1,las=y; for(int j=0;j<le;j++) g=mul(g,(y+n-x)%n,n),y=f(y,c,n); g=gcd(g,n); if(g==1) continue; if(g==n) { y=las,g=1; while(g==1) g=gcd((y+n-x)%n,n),y=f(y,c,n); } return g; } } } void Pollard_Rho(ll n,ll &a) { if(n<a) return; while(!(n&1)) n>>=1;a=2; if(n==1||Miller_Rabin(n)) {a=a>n?a:n;return;} ll d=Find(n); while(d==n) d=Find(n); if(d<n/d) d=n/d; Pollard_Rho(d,a),Pollard_Rho(n/d,a); }真实例题
CRT合并
考虑合并方程
\[ \left\{ \begin{aligned} x\equiv a_1\pmod {b_1}\\ x\equiv a_2\pmod {b_2} \end{aligned} \right. \]
由方程\(1\)，\(x=a_1+x_1b_1\)
代入方程\(2\)
\[ x_1b_1-x_2b_2=a_2-a_1 \]
等价于解同余方程
\[ ax+by=c\\ a=b_1,b=-b_2,c=a_2-a_1 \]
注意我们需要稍微调整使\(a,b,c\)非负
设\(d=\gcd(a,b)\)，同余方程有解当且仅当\(d|c\)，考虑解
\[ ax+by=d \]
假设解出一个特解为\(x_0,y_0\)，那么原方程有特解\(x_0'=x_0\frac{c}{d},y_0'=y_0\frac{c}{d}\)
考虑通解
\[ ax_0'+S+by_0'-S=c\\ a(x_0'+\frac{S}{a})+b(y_0'-\frac{S}{b})=c \]
那么正整数\(S\)最小取值为\(\operatorname{lcm}(a,b)\)
那么有通解\(x=x_0'+k\frac{\operatorname{lcm}(a,b)}{a}=x_0+k\frac{b}{\gcd(a,b)},k\in \mathbb Z\)
把\(a=b_1,b=-b_2,c=a_2-a_1\)代回去然后带到\(x=a_1+x_1b_1\)里面，可以得到
\[ \begin{aligned} x&=a_1+b_1(x_0'+k\frac{b_2}{\gcd(b_1,b_2)}),k\in \mathbb Z\\ &=a_1+b_1x_0'+k\operatorname{lcm}(b_1,b_2),k\in \mathbb Z \end{aligned} \]
这就是合并后的解是在\(\bmod \operatorname{lcm}(b_1,b_2)\)下的原因
要注意几个细节：
- 调整exgcd时的正负
- 从\(ax+by=d\)的特解到\(ax+by=c\)的特解的时候是对\(\frac{\operatorname{lcm}(a,b)}{a}\)取模，并且这个模一定要取，否则爆ll
- 乘法都写一写龟速乘，并不会改变复杂度。
真实例题
BSGS
BSGS
前置知识：欧拉定理：若正整数\(a,p\)互质，那么有\(a^{\varphi(p)}\equiv 1\pmod p\)
给定\(a,p,b\)，其中\((a,p)=1\)，求最小自然数解\(x\)
\[ a^x\equiv b\pmod p \]
根据欧拉定理，\(a^x\)最多只有\(\varphi(p)\)个，所以我们可以只考虑\(x\in [0,\varphi(p))\)的情况。首先特判几种情况：
- \(b=1\)时，\(x=0\)
- \(a=0\)时，若\(b\not=0\)，则无解，否则有解\(x=1\)
- \(b=0\)时，若\(a\not=0\)，则无解，否则有解\(x=1\)
注意这个特判是按顺序的，因为我们这里规定\(0^0=1\)
然后考虑对\(x\)进行分块，具体的，可以对\(x\)进行拆分
\[ \begin{aligned} &a^{ti-k}\equiv b\pmod p\\ &a^{ti}\equiv b\times a^k\pmod p\\ \end{aligned} \]
取\(t=\sqrt n\)，然后枚举\(k\in [0,t-1]\)，并把所有的\(ba^k\)插入\(\mathcal {Hash}\)表中。然后枚举\(i\in [1,t]\)，在\(\mathcal{Hash}\)表中查询\(a^{ti}\)，看是否有\(k\)可以匹配上。
exBSGS
给定\(a,p,b\)，求最小自然数解\(x\)
\[ a^x\equiv b\pmod p \]
一样，首先考虑特判几种情况：\(b=1\)时，\(x=0\)
\(a=0\)时，若\(b\not=0\)，则无解，否则有解\(x=1\)
\(b=0\)时
方程变成了这样
\[ a^x\equiv 0\pmod p \]
稍微转换一下
\[ k=\frac{a^x}{p},k\in \mathbb Z \]
发现直接枚举\(x\)就可以了，把\((a,p)\)除掉，然后看什么时候\(p=1\)，就返回\(x\)，如果\((a,p)=1\)，那么直接无解。
少了互质的条件，考虑构造互质。
和\(b=0\)时一个思路，考虑先枚举一部分\(x\)
\[ a\times a^{x-1}\equiv b\pmod p \]
稍微动一动式子
\[ a\times a^{x-1}-yp=b \]
根据裴蜀定理，如果\((a,p)\nmid b\)，那么返回无解否则把整个式子除\((a,p)\)，变成
\[ \frac{a}{(a,p)}a^{x-1}-y\frac{p}{(a,p)}=\frac{b}{(a,p)}\\ \frac{a}{(a,p)}a^{x-1}\equiv \frac{b}{(a,p)}\pmod {\frac{p}{(a,p)}} \]
考虑令\(A=a,B=\frac{b}{(a,p)},C=\frac{a}{(a,p)},P=\frac{p}{(a,p)}\)那么式子成了
\[ CA^x\equiv B\pmod P \]
继续递归即可，如果某一步\((A,P)=1\)，就用普通BSGS搞一下\(C\)显然不影响我们进行普通\(BSGS\)
注意一点，如果某一步\(C=B\)了，直接返回枚举的\(x\)，因为我们需要求最小自然数解。
枚举的\(x\)也就是递归深度是\(\log\)级别的，所以复杂度没问题
二次剩余
给定正整数\(a\)和奇素数\(p\)，求解同余方程
\[ x^2\equiv a\pmod p \]
欧拉准则：
若\(a^{\frac{p-1}{2}}\equiv 1\pmod p\)，那么\(a\)是\(\bmod p\)下的二次剩余
若\(a^{\frac{p-1}{2}}\equiv -1\pmod p\)，那么\(a\)不是\(\bmod p\)下的二次剩余
这里要判掉\(p|a\)的情况。
证明：
- 若是二次剩余，因为\((x,p)=1\)，根据费马小定理有\(x^{p-1}\equiv 1\pmod p\)，所以\(a^{\frac{p-1}{2}}\equiv 1\pmod p\)
- 若不是二次剩余，根据逆元的相互性，对同余方程\(ij\equiv x\pmod p\)，对\(i\in[1,p-1]\)，\(ij\)是一一相互配对的，那么就有\(a^{\frac{p-1}{2}}=(p-1)!\equiv -1 \pmod p\)，后面一步用了威尔逊定理。
判断有解后，考虑如何求出这个解。
设\(p-1=2^t\times s\)
若\(t=1\)，那么
\[ x\equiv \sqrt a=\sqrt{a\cdot a^{\frac{p-1}{2}}}=a^{\frac{s+1}{2}}\pmod p \]若\(t>1\)
设\(x_{t-1}=a^{\frac{s+1}{2}}\)，然后我们有
\[ (a^{-1}\cdot x^2_{t-1})^{2^{t-1}}=a^{\frac{p-1}{2}}\equiv 1\pmod p \]
然后我们把前面的东西拿出来，定义成一个类似与单位根的东西
\[ \epsilon_i=a^{-1}\cdot x_i^2 \]
可以发现，最后要求的是下式中的\(x_0\)
\[ (a^{-1}\cdot x^2_0)^{2^{0}}\equiv 1\pmod p \]
考虑\(x_{t-k}\)向\(x_{t-k-1}\)进行递推当\(k=1\)时，我们有\(x_{t-1}=a^{\frac{s+1}{2}},\epsilon_{t-1}=a^s\)
考虑从\(\epsilon_{t-k}\)入手向下推
我们知道当\(x^2\equiv 1\pmod p\)时，有\(x\equiv \pm 1\pmod p\)
所以有
\[ \epsilon_{t-k}^{2^{t-k-1}}\equiv \pm1\pmod p \]当\(\epsilon_{t-k}^{2^{t-k-1}}\equiv 1\pmod p\)时
我们发现直接让\(x_{t-k-1}=x_{t-k},\epsilon_{t-k-1}=\epsilon_{t-k}\)就可以了
当\(\epsilon_{t-k}^{2^{t-k-1}}\equiv -1\pmod p\)时
考虑找到一个整数\(\lambda\)，使得\(x_{t-k-1}=\lambda \cdot x_{t-k}\)
已知
\[ \epsilon_{t-k-1}^{2^{t-k-1}}=[a^{-1}\cdot(\lambda\cdot x_{t-k})^2]^{2^{t-k-1}}\equiv 1\pmod p \]
把里面的\(\lambda\)单独提出来，然后把\(\epsilon_{t-k}^{2^{t-k-1}}\equiv -1\pmod p\)带进去，可以得到
\[ (\lambda^2)^{2^{t-k-1}}=\lambda^{2^{t-k}}\equiv-1\pmod p \]
考虑方程\(x^2\equiv b\pmod p\)，若\(b\)在\(p\)意义下无二次剩余，那么根据欧拉准则有\(b^{\frac{p-1}{2}}\equiv -1 \pmod p\)而\(b\in[1,p-1]\)中，有无二次剩余的数的个数刚好相等（这个的证明思路大概是每个二次剩余都有不同的两个解），就是说这里面一半的数有二次剩余，一半没有，于是我们可以直接随机\(b\)找到它，期望次数是\(2\)，直接在最开始找一下就可以。
现在有了\(b\)，我们知道
\[ b^{\frac{p-1}{2}}=b^{2^{t-1}s}=b^{2^{t-k}2^{k-1}s}\equiv -1\pmod p \]
可以得到
\[ \lambda\equiv b^{2^{k-1}s}\pmod p \]
可以发现，复杂度大概是\(O(t^2)\)也就是\(O(\log^2 p)\)的
Timus Online Judge 1132. Square Root
upt:发现好像有一个\(\log\)的做法，自闭了
高次剩余
- 这个东西和二次剩余比起来简直小清新至极（但是我还是只学了奇素数版本的，主要是任意模数的听说写起来太长了）
- 参考资料：po姐的，原根与指标v.1.0.1，可以在U裙下载
考虑求解关于\(x\)的同余方程，其中\(P\)为奇素数（如无特殊说明，以下\(p,P\)全部为奇素数）
\[ x^A\equiv B\pmod P \]
两边取对数
\[ A\ln_? x_\equiv B\pmod ? \]
差不多是这个意思，先引入一些别的东西。
阶
定义关于\(x\)的同余方程\(a^x\equiv 1\pmod p\)的最小正整数解为\(a\)在\(\bmod p\)意义下的阶，记为\(\delta_p(a)\)
原根
若\(\delta_p(a)=\varphi(p)\)，那么称\(a\)为\(p\)的一个原根。
原根的一个重要性质
令\(g\)为\(p\)的一个原根，那么\(g^0,g^1,\dots,g^{\delta_p(g)-1}\)对\(\pmod p\)两两不同余
证明：考虑反证法，若存在\((j<i<\delta_p(g))\)满足
\[ g^i\equiv g^j\pmod p \]
那么有
\[ g^{i-j}\equiv 1\pmod p \]
则
\[ i-j<\delta_p(g) \]
与阶的定义矛盾。这个性质很有用，可以在后面作为离散对数的底。
原根的求法
设\(\varphi(p)=\prod d_i^{c_i}\)，根据定义，若\(g\)是原根，那么不存在
\[ g^{\frac{\varphi(p)}{d_i}}\equiv 1\pmod p \]
我们可以枚举原根，然后枚举\(d\)进行判断，因为原根一般都很小，所以求起来很快。
指标
若\(g\)为奇素数\(p\)的某个原根，若\(g^x\equiv a\pmod b\)，那么称\(x\)为以\(g\)为底\(\bmod m\)的一个指标，记做\(\operatorname {ind}_g(a)\)
这个定义和普通对数是相似的。
有了上面三个定义，我们回到之前的同余方程
\[ x^A\equiv B\pmod P \]
两边取对数
\[ A\operatorname {ind}_g(x)\equiv \operatorname{ind}_g(B)\pmod {\varphi (p)} \]
考虑使用BSGS解出\(\operatorname{ind}_g(B)\)，然后用扩欧求出\(\operatorname{ind}_g(A)\)的所有解，再代回去就可以了。
