数论ex
数学学得太差了补补知识点or复习
Miller-Rabin 和 Pollard Rho
-
Miller-Rabin
前置知识：
-
费马小定理
[a^{p-1}equiv 1pmod p,p is prime ] -
二次探测（mod奇素数下1的二次剩余）
[x^2equiv 1pmod pRightarrow x=1 or p-1 ]如果不是 (mod) 奇素数，二次剩余可能是更多的值
如果把费马小定理反过来用来检测一个数是否是素数，虽然是错的，但是反例比较少，如果配合二次探测进行测试并取多个a，可以把1e18内的所有数是否是质数判断出来。
具体的，取(a)为前(12)个质数((2,3,5,7,11,13,17,19,23,29,31,37))
然后用费马小定理进行测试，如果(a^{p-1}equiv 1pmod p)，就根据二次探测检验是否有(a^{(p-1)/2}equiv 1pmod 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; }这样做的复杂度是(12log^2 n)，其中因子2的个数是一个(log)，里面每次还做了一次快速幂，如果里面用了龟速乘就更完蛋了。
考虑从底向上做，就是说，先求出(x-1)除2到不能除的最底层的(x)，这样每次向上一层直接就是(x^2)了
考虑在某一层(x=1)，那么如果之前的一层的(x' ot=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-ymod n,n) ot=1,n)，那么(|x-y|)就是一个(n)的约数。
现在搞一个近似随机函数(f)，保证(f)在(mod 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(xymod n,n))
于是我们可以倍增(y)的步数，然后在每次倍增中，每(w)步（(w)可取(127,512)等值）取一次(gcd)，如果(gcd)值为(1)，说明这(gcd)步都没有值，继续向后找，如果(gcd)值为(n)，说明可能得到答案，就重新对这(w)求一遍，中间如果遇到(gcd ot=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合并
考虑合并方程
由方程(1)，(x=a_1+x_1b_1)
代入方程(2)
等价于解同余方程
注意我们需要稍微调整使(a,b,c)非负
设(d=gcd(a,b))，同余方程有解当且仅当(d|c)，考虑解
假设解出一个特解为(x_0,y_0)，那么原方程有特解(x_0'=x_0frac{c}{d},y_0'=y_0frac{c}{d})
考虑通解
那么正整数(S)最小取值为(operatorname{lcm}(a,b))
那么有通解(x=x_0'+kfrac{operatorname{lcm}(a,b)}{a}=x_0+kfrac{b}{gcd(a,b)},kin mathbb Z)
把(a=b_1,b=-b_2,c=a_2-a_1)代回去然后带到(x=a_1+x_1b_1)里面，可以得到
这就是合并后的解是在(mod 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 1pmod p)
给定(a,p,b)，其中((a,p)=1)，求最小自然数解(x)
[a^xequiv bpmod p ]根据欧拉定理，(a^x)最多只有(varphi(p))个，所以我们可以只考虑(xin [0,varphi(p)))的情况。
首先特判几种情况：
- (b=1)时，(x=0)
- (a=0)时，若(b ot=0)，则无解，否则有解(x=1)
- (b=0)时，若(a ot=0)，则无解，否则有解(x=1)
注意这个特判是按顺序的，因为我们这里规定(0^0=1)
然后考虑对(x)进行分块，具体的，可以对(x)进行拆分
取(t=sqrt n)，然后枚举(kin [0,t-1])，并把所有的(ba^k)插入(mathcal {Hash})表中。
然后枚举(iin [1,t])，在(mathcal{Hash})表中查询(a^{ti})，看是否有(k)可以匹配上。
-
exBSGS
给定(a,p,b)，求最小自然数解(x)
[a^xequiv bpmod p ]一样，首先考虑特判几种情况：
-
(b=1)时，(x=0)
-
(a=0)时，若(b ot=0)，则无解，否则有解(x=1)
-
(b=0)时
方程变成了这样
[a^xequiv 0pmod p ]稍微转换一下
[k=frac{a^x}{p},kin mathbb Z ]发现直接枚举(x)就可以了，把((a,p))除掉，然后看什么时候(p=1)，就返回(x)，如果((a,p)=1)，那么直接无解。
少了互质的条件，考虑构造互质。
和(b=0)时一个思路，考虑先枚举一部分(x)
[a imes a^{x-1}equiv bpmod p ]稍微动一动式子
[a imes a^{x-1}-yp=b ]根据裴蜀定理，如果((a,p) mid b)，那么返回无解
否则把整个式子除((a,p))，变成
[frac{a}{(a,p)}a^{x-1}-yfrac{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^xequiv Bpmod P ]继续递归即可，如果某一步((A,P)=1)，就用普通BSGS搞一下
(C)显然不影响我们进行普通(BSGS)
注意一点，如果某一步(C=B)了，直接返回枚举的(x)，因为我们需要求最小自然数解。
枚举的(x)也就是递归深度是(log)级别的，所以复杂度没问题
-
二次剩余
给定正整数(a)和奇素数(p)，求解同余方程
-
欧拉准则：
若(a^{frac{p-1}{2}}equiv 1pmod p)，那么(a)是(mod p)下的二次剩余
若(a^{frac{p-1}{2}}equiv -1pmod p)，那么(a)不是(mod p)下的二次剩余
这里要判掉(p|a)的情况。
证明：
- 若是二次剩余，因为((x,p)=1)，根据费马小定理有(x^{p-1}equiv 1pmod p)，所以(a^{frac{p-1}{2}}equiv 1pmod p)
- 若不是二次剩余，根据逆元的相互性，对同余方程(ijequiv xpmod p)，对(iin[1,p-1])，(ij)是一一相互配对的，那么就有(a^{frac{p-1}{2}}=(p-1)!equiv -1 pmod p)，后面一步用了威尔逊定理。
判断有解后，考虑如何求出这个解。
设(p-1=2^t imes s)
-
若(t=1)，那么
[xequiv sqrt a=sqrt{acdot 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 1pmod p ]然后我们把前面的东西拿出来，定义成一个类似与单位根的东西
[epsilon_i=a^{-1}cdot x_i^2 ]可以发现，最后要求的是下式中的(x_0)
[(a^{-1}cdot x^2_0)^{2^{0}}equiv 1pmod 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^2equiv 1pmod p)时，有(xequiv pm 1pmod p)
所以有
[epsilon_{t-k}^{2^{t-k-1}}equiv pm1pmod p ]-
当(epsilon_{t-k}^{2^{t-k-1}}equiv 1pmod p)时
我们发现直接让(x_{t-k-1}=x_{t-k},epsilon_{t-k-1}=epsilon_{t-k})就可以了
-
当(epsilon_{t-k}^{2^{t-k-1}}equiv -1pmod p)时
考虑找到一个整数(lambda)，使得(x_{t-k-1}=lambda cdot x_{t-k})
已知
[epsilon_{t-k-1}^{2^{t-k-1}}=[a^{-1}cdot(lambdacdot x_{t-k})^2]^{2^{t-k-1}}equiv 1pmod p ]把里面的(lambda)单独提出来，然后把(epsilon_{t-k}^{2^{t-k-1}}equiv -1pmod p)带进去，可以得到
[(lambda^2)^{2^{t-k-1}}=lambda^{2^{t-k}}equiv-1pmod p ]考虑方程(x^2equiv bpmod p)，若(b)在(p)意义下无二次剩余，那么根据欧拉准则有(b^{frac{p-1}{2}}equiv -1 pmod p)
而(bin[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 -1pmod p ]可以得到
[lambdaequiv 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^xequiv 1pmod p)的最小正整数解为(a)在(mod 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^iequiv g^jpmod p ]那么有
[g^{i-j}equiv 1pmod p ]则
[i-j<delta_p(g) ]与阶的定义矛盾。
这个性质很有用，可以在后面作为离散对数的底。
-
原根的求法
设(varphi(p)=prod d_i^{c_i})，根据定义，若(g)是原根，那么不存在
[g^{frac{varphi(p)}{d_i}}equiv 1pmod p ]我们可以枚举原根，然后枚举(d)进行判断，因为原根一般都很小，所以求起来很快。
-
-
指标
若(g)为奇素数(p)的某个原根，若(g^xequiv apmod b)，那么称(x)为以(g)为底(mod m)的一个指标，记做(operatorname {ind}_g(a))
这个定义和普通对数是相似的。
有了上面三个定义，我们回到之前的同余方程
两边取对数
考虑使用BSGS解出(operatorname{ind}_g(B))，然后用扩欧求出(operatorname{ind}_g(A))的所有解，再代回去就可以了。
单位根反演
- 老爹：要用魔法打败魔法
先丢一个式子
-
先复习一下复数的相关知识。
-
复数运算
[ ewcommand{i}{mathrm i} egin{aligned} &a+bi+c+di=(a+c)+(b+d)i\ &a+bi-(c+di)=(a-c)+(b-d)i\ &(a+bi)(c+di)=ac-bd+(ad+bc)i end{aligned} ]然后我们以实部为(x)轴，虚部为(y)轴把复数(z=a+bi)对于的坐标((a,b))放到平面直角坐标系中定义出复平面。
在复平面中的乘法，有模长相乘，幅角相加的法则，这个很好证，用三角函数算一下就可以了。
-
单位根（就是FFT里面那个）
我们发现，在单位圆上的复数相乘时具有优美的性质，模长始终为(1)。
定义
[z^n=1(n=1,2,3,dots) ]的复数根(z)为(n)次单位根，单位的(n)次根有(n)个
事实上，这(n)次单位根的(n)个根把单位元平均分成了(n)份，我们记第(k)份为(w_n^k)
有一些比较显然的运算性质：
[egin{aligned} &w_n^0=1\ &w_n^k=w_n^{kmod n}\ &(w_n^k)^i=w_n^{ki} end{aligned} ]
-
然后我们给出一个式子（如果学FFT认真的话，可以发现是用到这个式子了的）
这个式子很好证明，根据上面的运算性质，我们可以得到后面那个东西是等比数列的性质，然后你发现当(n|k)时所有项都是(1)，否则求一个等比数列发现分子是(0)
然后我们发现最开始的那个式子其实和这个式子差不多的，就是把最开始的单位根偏转了一个角度，对应的是指数上减去(l)就可以了，于是我们可以理解最开始的那个式子。
-
应用
我们知道，FFT做到最后(w)是没有虚部了的，所以可以直接用，正常情况下，我们比较难直接使用(w)
但是我们知道有个叫NTT的，用的是原根，代替了单位根。
很简单，我们只需要满足上面的运算性质就可以了。
我们知道(g^0,g^1,g^2,g^3,dots g^{varphi(p)-1} mod p)两两不同，且(g^0=1)
所以可以直接定义
[w_n^1equiv g^{frac{varphi(p)-1}{n}}pmod p ]这样一般题目给个(998244353)，就能直接用个什么(w_n^k)了
当然，这个推导还是要因题目而异的，不会只考一个单位根反演。
-
真实例题
类欧几里得 【simple】
- 式子推了三面，实在懒得抄上来，这里给出一些关键引理和最后结论。
- 这个东西思维不难，就是式子推起来有点恶心。
引理
-
[lfloorfrac{ai}{c} floor=lfloorfrac{amod c}{i} floor+ilfloorfrac{a}{c} floor ]
-
[X^2=-X+sum_{i=1}^Xi ]
-
[age bRightarrow a+1>b(a,bin mathbb Z) ]
-
[a>bRightarrow a>lfloor b floor(ain mathbb Z,bin mathbb R) ]
-
[lfloor a floor>bRightarrow a>b(ain mathbb R,bin mathbb Z) ]
都非常的显然把
结论
-
定义
[egin{aligned} &f(a,b,c,n)=sum_{i=0}^nlfloorfrac{ai+b}{c} floor\ &g(a.b,c,n)=sum_{i=0}^nlfloorfrac{ai+b}{c} floor^2\ &h(a.b,c,n)=sum_{i=0}^nilfloorfrac{ai+b}{c} floor end{aligned} ] -
结论
[egin{aligned} &f(a,b,c,n)=left{egin{aligned} &(n+1)lfloorfrac{b}{c} floor,a=0\ &f(amod c,b mod c,c,n)+frac{n(n+1)}{2}lfloorfrac{a}{c} floor+(n+1)lfloorfrac{b}{c} floor,age c or bge c\ &nlfloorfrac{an+b}{c} floor-f(c,c-b-1,a,lfloorfrac{an+b}{c} floor-1),a<c and b<c end{aligned} ight.\ &g(a,b,c,n)=left{ egin{aligned} &(n+1)lfloorfrac{b}{c} floor^2,a=0\ &g(amod c,b mod c,c,n)+frac{(2n+1)(n+1)n}{6}lfloorfrac{a}{c} floor^2+(n+1)lfloorfrac{b}{c} floor^2+2lfloorfrac{a}{c} floor h(amod c,b mod c,c,n)\&+2lfloorfrac{b}{c} floor f(amod c,b mod c,c,n)+n(n+1)lfloorfrac{a}{c} floor lfloorfrac{b}{c} floor,age c or bge c\ &n(lfloorfrac{an+b}{c} floor)(lfloorfrac{an+b}{c} floor+1)-2h(c,c-b-1,a,lfloorfrac{an+b}{c} floor-1)-2f(c,c-b-1,a,lfloorfrac{an+b}{c} floor-1)\&-f(a,b,c,n),a<c and b<c\ end{aligned} ight.\ &h(a,b,c,n)=left{egin{aligned} &frac{n(n+1)}{2}lfloorfrac{b}{c} floor,a=0\ &h(amod c,b mod c,n)+frac{(2n+1)(n+1)n}{6}lfloorfrac{a}{c} floor+frac{n(n+1)}{2}lfloorfrac{b}{c} floor,age c or bge c\ &frac{1}{2}[n(n+1)lfloorfrac{an+b}{c} floor-g(c,c-b-1,a,lfloorfrac{an+b}{c} floor-1)-f(c,c-b-1,a,lfloorfrac{an+b}{c} floor)]\ end{aligned} ight. end{aligned} ] -
这些推起来都不难，别被劝退了，就是麻烦了点。