数论算法 剩余系相关 学习笔记 (基础回顾,(ex)CRT,(ex)lucas,(ex)BSGS,原根与指标入门,高次剩余,Miller_Rabin+Pollard_Rho)
注:转载本文须标明出处。
原文链接https://www.cnblogs.com/zhouzhendong/p/Number-theory.html
数论算法 剩余系相关 学习笔记 (基础回顾,(ex)CRT,(ex)lucas,(ex)BSGS,原根与指标入门,高次剩余,Miller_Robin+Pollard_Rho)
本文概要
1. 基础回顾
2. 中国剩余定理 (CRT) 及其扩展
3. 卢卡斯定理 (lucas) 及其扩展
4. 大步小步算法 (BSGS) 及其扩展
5. 原根与指标入门
6. 从简到难解决高次剩余问题
7. 大素数判定与大数分解
8. 相关推荐
9. 鸣谢
基础回顾
欧几里得算法与扩展欧几里得算法
1 2 3 | int gcd( int x, int y){ return y?gcd(y,x%y):x; } |
1 2 3 4 5 6 7 8 9 | int exgcd( int a,ini b, int &x, int &y){ if (!b){ x=1,y=0; return a; } int res=exgcd(b,a%b,y,x); y-=(a/b)*x; return res; } |
快速幂
1 2 3 4 5 6 7 | int Pow( int x, int y, int mod){ int ans=1; for (;y;y>>=1,x=1LL*x*x%mod) if (y&1) ans=1LL*ans*x%mod; return ans; } |
欧拉函数
1 2 3 4 5 6 7 8 9 10 11 12 | int phi( int n){ int res=n; for ( int i=2;i*i<=n;i++) if (n%i==0){ res=res/i*(i-1); while (n%i==0) n/=i; } if (n>1) res=res/n*(n-1); return res; } |
欧拉定理与费马小定理
欧拉定理 :如果 $a,p$ 互质,则 $a^{\varphi(p)} \equiv 1 \pmod p$ 。
费马小定理 :当 $p$ 为质数的时候,则 $a^{p-1} \equiv 1 \pmod p$ 。
逆元
给定 $a,p$ ,求 $a^{-1} \pmod p$ 。
求法:如果 $a,p$ 不互质,显然无解。否则可以通过 扩展欧几里得 或 费马小定理/欧拉定理 来求。
当 $p$ 为质数的时候,运用费马小定理可以得到:
$$a^{-1}\equiv a^{p-2} \pmod p$$
中国剩余定理
问题模型
给定同余方程组
$$\begin{cases}x&\equiv&x_1&\pmod {p_1}\\x&\equiv&x_2&\pmod {p_2}\\ &&\vdots\\x&\equiv&x_n&\pmod {p_n}\end{cases}$$
求解 $x$ 在对 ${\rm lcm}(p_1,p_2,\cdots ,p_n)$ 取模的意义下的值,或判断无解。
模数两两互质的情况 (CRT)
我们考虑一种构造方法。
首先,设 $P=p_1p_2\cdots p_n={\rm lcm}(p_1,p_2,\cdots ,p_n)$ ,$A_i=\cfrac{P}{p_i}$ ,$B_i\equiv A_i^{-1} \pmod {p_i}$ 。
则,
$$x\equiv \sum_{i=1}^{n} x_iA_iB_i \pmod P$$
这就是著名的孙子定理,即中国剩余定理。
接下来我们来证明这个构造的正确性。
首先,由于 $p_i$ 与 $A_i$ 互质,所以必然存在逆元 $B_i$ 。
然后,由于 $p_j|A_i (j\neq i)$ ,所以 $(\sum_{i=1}^{n}x_iA_iB_i)\equiv x_jA_jB_j \pmod {p_j}$ 。
又由于 $B_j\equiv A_j^{-1} \pmod {p_j}$ ,所以 $x_jA_jB_j\equiv x_j \pmod {p_j}$ 。
所以构造成立。
不保证模数互质的情况 (exCRT)
接下来我们来解决扩展中国剩余定理。(下述证明摘自博主另一篇博文:https://www.cnblogs.com/zhouzhendong/p/exCRT.html)
我们只需要知道如何合并两个方程,就可以将原方程组逐个合并了,所以下面只介绍合并两个方程的方法。
$$\begin{eqnarray*}x&\equiv& a_1&\pmod {p_1}\\x&\equiv& a_2&\pmod {p_2}\end{eqnarray*}$$
则 $x=a_1+k_1p_1=a_2+k_2p_2$ 。
令 $p={\rm lcm}(p_1,p_2)$,需要求形同 $x\equiv a\pmod p$的一个式子。
$$\begin{eqnarray*}a_1+k_1p_1&\equiv&a_2+k_2p_2 &\pmod {p}\\a_2-a_1&\equiv&k_1p_1-k_2p_2 & \pmod {p}\end{eqnarray*}$$
我们先用 exgcd 求解方程 $p_1x+p_2y=\gcd(p_1,p_2)$。设 $g=\gcd(p_1,p_2),z=a_2-a_1$。
由于 $g|p$ ,所以,如果 $g|z$ 不成立,那么显然 $g+kp|z$ 也不成立,故无解。
否则有解。记 $t=\frac zg$ ,则 $k_1=tx_1,k_2=-tx_2$,于是,$a\equiv a_1+k_1p_1\equiv a_2+k_2p_2\pmod p$ ,单次合并完成。
扩展中国剩余定理代码
1 2 3 4 5 6 7 8 9 10 | bool CRT(LL w1,LL p1,LL w2,LL p2,LL &w,LL &p){ LL x,y,z=w2-w1,g=ex_gcd(p1,p2,x,y); if (z%g) return 0; LL t=z/g; x=Mul(x,t,p2/g); p=p1/g*p2; w=((w1+Mul(x,p1,p))%p+p)%p; return 1; } |
例题传送门
POJ2891 - 模板题
BZOJ2010
NOI2018Day2T1 - 经典题
51Nod1362
卢卡斯定理
定义
设 $n,m$ 为正整数, $p$ 为素数,若
$$n=a_kp^k+a_{k-1}p^{k-1}+\cdots a_1p+a_0$$
$$m=b_kp^k+b_{k-1}p^{k-1}+\cdots b_1p+b_0$$
所有的 $a_i,b_i$ 都为非负整数,则:
$$\binom{n}{m} \equiv \prod _{i=0}^{k} \binom {a_i}{b_i}\pmod p$$
或者用递归表示:
$$\binom{n}{m} \equiv \binom {\left\lfloor\frac{n}{p}\right\rfloor}{\left\lfloor\frac{m}{p}\right\rfloor}\binom{n \bmod p}{m \bmod p}$$
证明
下面首先摘录一段“冯志刚版《初等数论》第37页”的内容来证明。貌似也可以从组合意义方向yy。
由 $p$ 为素数,可知对 $1\leq j\leq p-1 $ ,都有 $\binom{p}{j}=\cfrac{p}{j}\binom{p-1}{j-1}\equiv 0 \pmod p $,于是,
$$\begin{eqnarray*}(1+x)^p&=&1+\binom{p}{1}x+\cdots + \binom{p}{p-1}x^{p-1}+x^p \\ &\equiv & 1 + x^p \pmod p \end{eqnarray*}$$
利用上述结果,可知
$$\begin{eqnarray*}(1+x)^a&=&(1+x)^{a_0}\times((1+x)^p)^{a_1}\times \cdots \times ((1+x)^{p^k})^{a_k} \\ &\equiv& (1+x)^{a_0} \times (1+x^p)^{a_1} \times \cdots \times (1+x^{pk})^{a_k} \pmod p \end{eqnarray*}$$
对比两边 $x^m$ 项的系数(用二项式定理及 $p$ 进制数的性质)可得
$$\binom{n}{m} \equiv \prod _{i=0}^{k} \binom {a_i}{b_i}\pmod p$$
代码
1 2 3 4 5 | LL Lucas(LL n,LL m){ if (n<mod&&m<mod) return C(n,m); return C(n%mod,m%mod)*Lucas(n/mod,m/mod); } |
扩展卢卡斯
考虑模数不是质数的情况。
首先对模数分解质因数,令 $P=\prod p_i^{k_i}$ ,其中 $p_i$ 都是指数。我们只需要对于每一个 $p_i^{k_1}$ 为模数的时候求出答案,最后 CRT 合并就可以了。
接下来考虑对于 $p^k$ 形式的模数取模的做法。其中 $p$ 为质数。
由于 $ap^x$ 形式的数在 对于 $p^k$ 取模意义下没有逆元,所以我们要对于这一类数特殊处理。扩展卢卡斯的关键在于除去和统计阶乘中含 $p$ 因子的个数。
首先,我们考虑计算除去含 $p$ 因子后的阶乘在对于 $p^k$ 取模意义下的值。
考虑一个数列 $1$~$n$ :
$$1,2,\cdots p-1,p,p+1,\cdots ,2p,\cdots ,kp, \cdots ,\left\lfloor\frac{n}{p} \right \rfloor p, \cdots ,n$$
其中,所有不为 $p$ 倍数的数都直接对于我们要算的东西有直接贡献,我们直接把它累乘到答案中;一般来说 $p^k$ 都不会很大,所以我们可以枚举 $1$ 到 $p^k$ 中与 $p$ 互质的数累乘,超过 $p^k$ 之后,就相当于从 $0$ 开始循环了。
对于 $p$ 的倍数的数,我们将它们都除掉一个 $p$ ,然后递归调用上述过程,直到没有这种数了为止。时间复杂度 $O(p^k)$ 。
然后,我们统计组合数中包含的 $p$ 因子的个数。这个比较容易。
考虑 $1$~$n$ 的自然数中,至少含有一个因子 $p$ 的数的个数为 $\left\lfloor \frac {n}{p} \right\rfloor $ 。然后除去这些数的因子 $p$ ,不断求 $1$~$\left\lfloor \frac {n}{p} \right\rfloor $ 之间的至少含有一个因子 $p$ 的数的个数即可。
代码
代码貌似还不短啊。但是比起后面的高次剩余算不了什么了。
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 | LL Pow(LL x,LL y,LL mod){ if (y==0) return 1LL; LL xx=Pow(x,y/2,mod); xx=xx*xx%mod; if (y&1LL) xx=xx*x%mod; return xx; } void ex_gcd(LL a,LL b,LL &x,LL &y){ if (!b) x=1,y=0; else ex_gcd(b,a%b,y,x),y-=a/b*x; } LL Inv(LL X,LL mod){ if (!X) return 0; LL a=X,b=mod,x,y; ex_gcd(a,b,x,y); x=(x%b+b)%b; return x; } LL ex_lucas(LL n,LL pi,LL pk){ if (!n) return 1LL; LL ans=1; for (LL i=2;i<=pk;i++) if (i%pi) ans=ans*i%pk; ans=Pow(ans,n/pk,pk); for (LL i=2;i<=n%pk;i++) if (i%pi) ans=ans*i%pk; return ans*ex_lucas(n/pi,pi,pk)%pk; } LL C(LL n,LL m,LL pi,LL pk){ if (m>n) return 0; LL a=ex_lucas(n,pi,pk),b=ex_lucas(m,pi,pk),c=ex_lucas(n-m,pi,pk); LL k=0,ans; for (LL i=n;i;i/=pi,k+=i); for (LL i=m;i;i/=pi,k-=i); for (LL i=n-m;i;i/=pi,k-=i); ans=a*Inv(b,pk)%pk*Inv(c,pk)%pk*Pow(pi,k,pk)%pk; return ans*(P/pk)%P*Inv(P/pk,pk)%P; } LL C(LL n,LL m){ LL ans=0; for ( int i=1;i<=cnt;i++) ans=(ans+C(n,m,px[i],Pow(px[i],py[i],P+1)))%P; return ans; } |
例题
BZOJ2142
大步小步算法 (BSGS)
问题模型
求解 $x$ ,满足 $A^x\equiv B \pmod p$ 。$1\leq A,B,p\leq 10^9$
BSGS
BSGS 算法可以解决的是 $A,p$ 互质的情况。
首先确定一个阀值 $M$ (一般取 $\sqrt p$ 的近似值),则最终的解都可以表示成 $aM-b$ 的形式。令
$$A^{aM-b}\equiv B \pmod p$$
则可以得到
$$A^{aM}\equiv B\times A^b \pmod p$$
于是,我们只需要先枚举 $b$ 的值,把对应的 $B\times A^b$ 扔到一个 $Map$ 里面,然后再枚举左侧的 $a$ ,查询满足上式条件的 $b$ 的值即可。
时间复杂度 $O(\sqrt p \log p)$ 。用无序Map (unordered_map) 可以省掉 $\log$ ,但是常数较大。推荐的方法是手写 Hash 表,时间复杂度 $O(\sqrt p)$ 。
代码
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 | struct hash_map{ static const int Ti=233,mod=1<<16; int cnt,k[mod+1],v[mod+1],nxt[mod+1],fst[mod+1]; int Hash( int x){ int v=x&(mod-1); return v==0?mod:v; } void clear(){ cnt=0; memset (fst,0, sizeof fst); } void update( int x, int a){ int y=Hash(x); for ( int p=fst[y];p;p=nxt[p]) if (k[p]==x){ v[p]=a; return ; } k[++cnt]=x,nxt[cnt]=fst[y],fst[y]=cnt,v[cnt]=a; return ; } int find( int x){ int y=Hash(x); for ( int p=fst[y];p;p=nxt[p]) if (k[p]==x) return v[p]; return 0; } int &operator [] ( int x){ int y=Hash(x); for ( int p=fst[y];p;p=nxt[p]) if (k[p]==x) return v[p]; k[++cnt]=x,nxt[cnt]=fst[y],fst[y]=cnt; return v[cnt]=0; } }Map; int BSGS( int A, int B, int P){ int M=max(( int )(0.8* sqrt (1.0*P)),1),AM=Pow(A,M,P); Map.clear(); for ( int b=0,pw=B;b<M;b++,pw=1LL*pw*A%P) Map.update(pw,b+1); for ( int a=M,pw=AM;a-M<P;a+=M,pw=1LL*pw*AM%P){ int v=Map.find(pw); if (v) return a-(v-1); } return -1; } |
把 Hash 表改成 unordered_map 就只有这么一点点代码了:
1 2 3 4 5 6 7 8 9 10 11 | unordered_map < int , int > Map; int BSGS( int A, int B, int P){ int M=max(( int )(0.8* sqrt (1.0*P)),1),AM=Pow(A,M,P); Map.clear(); for ( int b=0,pw=B;b<M;b++,pw=1LL*pw*A%P) Map[pw]=b+1; for ( int a=M,pw=AM;a-M<P;a+=M,pw=1LL*pw*AM%P) if (Map[pw]) return a-(Map[pw]-1); return -1; } |
扩展BSGS
不保证 $A,p$ 互质。
令 $g=\gcd(A,p)$ ,则:
$$\frac{A}{g} A^{x-1} \equiv \frac {B}{g} \pmod {\frac{p}{g}}$$
当 $B$ 不是 $g$ 的倍数的时候,显然无解。
于是又得到了一个新的问题,性质与原问题一样。于是可以不断递归求解,直到 $g=1$ 为止。
最终会得到一个形如
$$\frac {A^k} {d} A^{x-k} \equiv \frac{B}{d} \pmod {\frac {p}{d}}$$
的方程。
其中 $d$ 为各个 $g$ 的乘积。
于是问题转化成了求解方程
$$ a A^x \equiv b \pmod {p^\prime}$$
于是 BSGS 就可以解决的问题了,最终把 $k$ 加回去就可以了。
注意,上述做法不能解决 $x<k$ 的情况,对于这种情况,我们只需要特殊处理一下就可以了。
代码
此处略去 Hash 表。
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 | int ExBSGS( int A, int B, int P){ A%=P,B%=P; int k=0,v=1; while (1){ int g=gcd(A,P); if (g==1) break ; if (B%g) return -1; k++,B/=g,P/=g,v=1LL*v*(A/g)%P; if (v==B) return k; } if (P==1) return k; int M=max(( int ) sqrt (1.0*P),1),AM=Pow(A,M,P); Map.clear(); for ( int b=0,pw=B;b<M;b+=1,pw=1LL*pw*A%P) Map.update(pw,b+1); for ( int a=M,pw=1LL*v*AM%P;a-M<P;a+=M,pw=1LL*pw*AM%P){ int v=Map.find(pw); if (v) return a-(v-1)+k; } return -1; } |
例题
BZOJ2480 - 扩展bsgs模板题
原根与指标入门
声明:这里摘取了部分 PoPoQQQ 的课件《原根与指标 v1.0.1 .ppt》的内容,感谢 PoPoQQQ。读者如果想要看这个课件,请自行移步 UOJ 群文件中下载。
阶
定义
对于一个正整数 $a$ ,求满足 $a^x \equiv 1 \pmod m$ 的最小正整数 $x$ 。保证 $a,m$ 互质。这个最小正整数 $x$ 称作 $a$ 在对 $m$ 取模意义下的阶,记做 ${\rm ord}_m(a)$ ,在 $m$ 的值十分明确的时候,可以记做 ${\rm ord}(a)$ 。
性质
1. $a^n\equiv 1 \pmod m$ 的充要条件是 ${\rm ord}_m(a) | n$ 。推论: ${\rm ord}_m(a) | \varphi (m) $ 。
2. 若 $a\equiv b \pmod m$ ,则 ${\rm ord} _m(a) = {\rm ord}_m(b)$ 。
3. 若 $a^n\equiv a^i \pmod m $,则 $n\equiv i \pmod {{\rm ord}_m(a)}$ 。
4. 令 $n= {\rm ord} _m (a) $ ,则 $a^0,a^1,\cdots ,a^{n-1}$ 对 $m$ 取模两两不同。
由于上述性质的证明比较简单,所以这里就不写证明了。
原根
定义
若 $g$ 是模 $p$ 意义下的原根,则 $g$ 满足 ${\rm ord}_p(g) = \varphi(p)$ 。
性质
1. 模 $p$ 意义下存在原根,当且仅当 $p$ 是如下形式的数: $2,4,x^a,2x^a$ 。($x$ 为奇素数, $a$ 为正整数)
2. 当 $p$ 为奇素数时,模 $p$ 意义下的原根个数为 $\varphi(\varphi(p))$ 。
3. 对于质数 $p$ , $\varphi (p)=p-1$ ,将 $p-1$ 分解质因数,得到 $p-1=\prod {p_i^{q_i}}$ ,则正整数 $g$ 是模 $p$ 意义下的原根的充分必要条件是:对于所有 $i$ ,$g^{\frac{p-1}{p_i}}\not\equiv 1 \pmod p$ 。 证明: 充分性很显然。必要性:首先考虑阶的第 1 点性质,可以得知 ${\rm ord}_p(g)| p-1$ ,那么,如果这个值比 $p-1$ 小,必然可以找到一个 $i$ ,使得 ${\rm ord}_p(g) | \frac{p-1}{p_i} $ ,那么 $g^{\frac{p-1}{p_i}}\equiv 1 \pmod p$ ,故 $g$ 不是原根,否则,说明 ${\rm ord }_p(m) = p-1 =\varphi(p)$ ,$g$ 是原根。
如何求原根
由于原根一般都很小,所以我们直接暴搜原根!
具体多小?根据 PoPoQQQ 的课件,在 $10^9$ 范围内,最小的原根大小大约不超过 $n^{0.25}$ 。
代码
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | bool Get_g_Check( int P, int x){ for ( int i=1;i<=Fac_tot;i++) if (Pow(x,(P-1)/Fac_p[i],P)==1) return 0; return 1; } int Get_g( int P){ Fac_tot=0; int v=P-1; for ( int i=1;prime[i]*prime[i]<=v&&i<=pcnt;i++) if (v%prime[i]==0){ Fac_p[++Fac_tot]=prime[i]; while (v%prime[i]==0) v/=prime[i]; } if (v>1) Fac_p[++Fac_tot]=v; for ( int i=2;;i++) if (Get_g_Check(P,i)) return i; return -1; } |
指标
定义
设 $g$ 是模 $p$ 意义下的一个原根,若 $a,p$ 互质,则存在一个唯一的正整数 $x$ 满足 $g^x\equiv a \pmod m $ 。这个 $x$ 就叫做 “以 $g$ 为底的 $a$ 对模的一个指标",记作 $x={\rm ind}_g(a) $ 。
求法
直接 BSGS 。
高次剩余
问题模型
求解方程
$$x^A \equiv B \pmod p$$
$1\leq A,B,p \leq 10^9$
当 $p$ 为质数时
传送门 - 51Nod1038题意与代码
首先,我们可以找到模 $p$ 意义下的原根 $g$ 。
令 $t={\rm ind}_g(B)$ ,设 $x=g^y$ ,则
$$g^{Ay} \equiv g^t \pmod p$$
即
$$Ay \equiv t \pmod {\varphi (p)}$$
把 $y$ 解出来即可得到 $x$ 。
代码见传送门。
当 $p$ 为奇数时
传送门 - BZOJ2219题意与代码
注:这里这道 bzoj 题只要求输出解的个数。
考虑对 $p$ 分解质因数。
假设 $p=\prod p_i^{q_i}$ ,$p_i$ 为奇素数。那么我们显然可以对于每一个 $p_i$ 分开求解,最后可以通过中国剩余定理得到:在原模数意义下的解的总数,为所有 模 $p_i^{q_i}$ 意义下的解的个数 的乘积。
考虑求解形如 $x^A \equiv B \pmod {p^k}$ 的方程。
如果 $\gcd(B,p^k) = p^k$
那么,只需要保证 $p^{\left\lceil\frac kA\right\rceil}|x$ 即可。
令 $base = p^{\left\lceil\frac kA\right\rceil}$ ,则解的形式为 $\alpha\cdot base$ ,解的个数为 $p^{k-\left\lceil\frac kA\right\rceil} $ 。
如果 $\gcd(B,p^k) = p^q \ \ \ (0\leq q<k)$
则显然有 $A | q$ 。
令 $B^\prime = \frac B {p^q}$ ,则可以把原方程转化成求
$$p^q y^A \equiv p^q B^\prime \pmod {p^k}$$
即
$$y^A \equiv B^\prime \pmod {p^{k-q}}$$
由于 $\gcd(B^\prime , p^{k-q})=1$ ,那么求解 $y$ 显然可以用之前介绍的方法来得到。
最终, $x \equiv y \times p^{\frac qA} \pmod {p^k}$ 。于是,如果 $y$ 在模 $p^{k-q}$ 意义下的解的个数为 $ans$ ,那么 $x$ 在模 $p^k$ 意义下的解的个数为 $p^{q-\frac qA}\times ans$ 。
代码见传送门。
当 $p$ 没有限制时
传送门 - 51Nod1123题意与代码
其他都一样,在求解 $x^A \equiv B \pmod {p^k}$ 这一类问题时,当 $p=2$ 时,问题会比较特殊。
考虑到 $x^A \equiv B \pmod {2^k}$ 的充要条件是 $x^A \equiv B \pmod {2^{k-1}}$ 。
反过来,如果 $x^A \equiv B \pmod {2^{i-1}}$ ,那么只可能导致 $x^A \equiv B \pmod {2^i}$ 或者 $(x+2^{i-1})^A \equiv B \pmod {2^i}$ 。只需要写个 dfs 就好了。
但是由于这题要求出具体的解,所以代码有点长。
代码见传送门。
小结
高次剩余的求解中,要用到之前提到过的几乎所有算法。数论算法在我们眼中突然变成了码农题,240行高次剩余,瞬间虐爆数据结构。
”告辞“! 剩余!
如果您还想来一道模板题的话
51Nod1039
大素数判定和大整数分解算法 - Miller_Rabin & 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 | #include <bits/stdc++.h> using namespace std; typedef long long LL; typedef unsigned long long ULL; typedef long double LD; namespace Pollard_Rho{ int prime[9]={2,3,5,7,11,13,17,19,23}; ULL RR; int Pcnt; LL p[70]; vector <LL> res; LL R(LL mod){ return (RR+=4179340454199820289LL)%mod; } LL Mul(LL x,LL y,LL mod){ LL d=(LL) floor ((LD)x*y/mod+0.5); LL res=x*y-d*mod; if (res<0) res+=mod; return res; } LL Pow(LL x,LL y,LL mod){ LL ans=1%mod; for (;y;y>>=1,x=Mul(x,x,mod)) if (y&1) ans=Mul(ans,x,mod); return ans; } bool Miller_Rabin(LL n){ if (n<=1) return 0; for ( int i=0;i<9;i++) if (n==prime[i]) return 1; LL d=n-1; int tmp=0; while (!(d&1)) d>>=1,tmp++; for ( int i=0;i<9;i++){ LL x=Pow(prime[i],d,n),p=x; for ( int j=1;j<=tmp;j++){ x=Mul(x,x,n); if (x==1&&p!=1&&p!=n-1) return 0; p=x; } if (x!=1) return 0; } return 1; } LL f(LL x,LL c,LL mod){ return (Mul(x,x,mod)+c)%mod; } LL gcd(LL x,LL y){ return y?gcd(y,x%y):x; } LL Get_Factor(LL c,LL n){ LL x=R(n),y=f(x,c,n),p=n; while (x!=y&&(p==n||p==1)){ p=gcd(n,max(x-y,y-x)); x=f(x,c,n); y=f(f(y,c,n),c,n); } return p; } void Pollard_Rho(LL n){ if (n<=1) return ; if (Miller_Rabin(n)){ res.push_back(n); return ; } while (1){ LL v=Get_Factor(R(n-1)+1,n); if (v!=n&&v!=1){ Pollard_Rho(v); Pollard_Rho(n/v); return ; } } } void work(LL n){ res.clear(); Pollard_Rho(n); } } |
相关推荐
基础数论函数相关 - 莫比乌斯反演 - 杜教筛 - https://www.cnblogs.com/zhouzhendong/p/Number-theory-Function.html
快速傅里叶变换/快速数论变换 入门与经典例题 - https://www.cnblogs.com/zhouzhendong/p/Fast-Fourier-Transform.html
鸣谢
冯志刚版《初等数论》
yyb - BSGS - http://www.cnblogs.com/cjyyb/p/8810050.html
PoPoQQQ 原根与指标 v1.0.1 .ppt
感谢 51Nod 和 BZOJ 提供大量例题来源。
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 为什么说在企业级应用开发中,后端往往是效率杀手?
· 用 C# 插值字符串处理器写一个 sscanf
· Java 中堆内存和栈内存上的数据分布和特点
· 开发中对象命名的一点思考
· .NET Core内存结构体系(Windows环境)底层原理浅谈
· 为什么说在企业级应用开发中,后端往往是效率杀手?
· 本地部署DeepSeek后,没有好看的交互界面怎么行!
· DeepSeek 解答了困扰我五年的技术问题。时代确实变了!
· 趁着过年的时候手搓了一个低代码框架
· 推荐一个DeepSeek 大模型的免费 API 项目!兼容OpenAI接口!