标签：return 因数分解 ll register Pollard Rho 浅谈 mod
浅谈质因数分解
->part 1：
算数基本定理：
任何一个大于1的正整数都能唯一分解为有限个质数的乘积，可写作：
\[N=\prod_{i=1}^m p_i^ {c_i}\]
其中\(c_i\)都是正整数，\(p_i\)都是质数，且满足\(p_1<p_2<…<p_m\)
->part 2：
分解方法：
- 试除法
结合质数判定的“试除法”和质数筛选的“\(Eratosthenes\) 筛法”,我们可以扫描 \(2-\sqrt N\)的每个数\(x∈Z\)，若\(x\)能整除\(N\)，则从\(N\)中除掉所有质因子\(x\)，同时累计除去\(x\)的个数。
上述算法保证一个合数的因子一定在扫描到这个合数之前除掉了，所以这个过程中所得到的能整除\(N\)的一定是质数。时间复杂度为\(O(\sqrt N)\)。
值得注意的一点是，若\(N\)没有被任何\(2-\sqrt N\)的整数整除，则\(N\)是质数，无需分解。
- Pollard-Rho 算法
试除法的时间复杂度在\(O(\sqrt N)\)而 \(Pollard Rho\) 算法的理论实践复杂度在\(O(\sqrt[4]{N})\)，这已经大大提高了程序效率。
学习 \(Pollard Rho\) 算法前，需要掌握几个知识点：
1. Miller-Rabin 算法
- 学习要求：掌握费马小定理、二次探测定理
因此，要对同余相关知识有所了解，详见同余相关及中国剩余定理
- 费马小定理：
若\(p\)为质数，则
\[ ∀a∈Z,a^p \equiv a\pmod{p}\]
费马小定理的证明可以参考《算法竞赛进阶指南》，读者也可以自己尝试证明或者百度。
上述定理对于同个\(p\)、多个底数 \(a\)，若不满足上式，则\(p\)为合数。
- 二次探测定理：
若\(p\)为质数，
且 \(x^2 \equiv 1\pmod{p}\)， 则\(x \equiv 1 \pmod{p}\) 或 \(x \equiv p-1 \pmod{p}\)
二次探测定理就用作提高程序准确性：
我们对于一个数\(p\)且满足：$ a^{p-1} \equiv 1\pmod{p}$[费马小定理]
那么对于\(2|p-1\)，可以得到\((x^\frac{p-1}{2})^2 \equiv 1 \pmod{p}\)，我们在费马小定理的基础上讨论$x^\frac{p-1}{2} \mod{p} \(的值(以下记为\)K$)即可：
- 若\(K ≠ 1\) 且 $K ≠p-1 \(，则\)p$是合数 \(return \ false\);
- 若\(K=1\)，则继续递归\(x^\frac{p-1}{2}\);
- 若\(K=p-1\),不能利用二次探测定理继续递归，说明目前无法验证\(p\)是合数，\(return \ true\);
(以上部分为引用内容)
这样我们多取几个\(x\)就能够提高正确率。
以上内容是为了实现Miller-Rabin 算法，这是一个高效\((O(logN))\)的判断素数的随机算法。
昨天两个毒瘤数据还是卡了我很久，感谢 @雪中舞 巨佬的帮助
然后拿两个质数 \(61\) 和 \(24251\) 以及三个随机数去\(check\)，这样就基本不会被卡了。
然而我三个随机数还是被卡了，因为某种意外，详见下面代码；于是我写了四个随机数，这样就稳稳AC了。
//Miller_Rabin 算法
inline ll qmul(ll x, ll y, ll mod) { //快速乘
return (x*y-(ll)((long double)x/mod*y)*mod+mod)%mod;
}
inline ll qpow(ll x, ll m, ll mod) { //快速幂
register ll ans = 1;
while(m) {
if(m & 1)
ans = qmul(ans, x, mod) % mod;
x = qmul(x, x, mod) % mod;
m >>= 1;
}
return ans;
}
inline bool FMLT(ll x, ll p) { //费马小定理
return qpow(x, p - 1, p) == 1;
}
inline bool Miller_Rabin(ll x, ll p) {
if(!FMLT(x, p)) return 0;
register ll k = p - 1;
for(;!(k & 1);) {
k >>= 1;
ll K = qpow(x, k, p);
if(K != 1 && K != p - 1) return 0;
if(K == p - 1) return 1;
}
return 1;
}
inline bool check(ll p) { // 判断质数
if(p == 1) return 0;
ll t1 = rand(), t2 = rand(), t3 = rand(), t4 = rand();
if(p == 61 || p == 24251 || p == t1 || p == t2 || p == t3 || p == t4) return 1;
/*这行写成
if(p == 2 || p == 3 || p == t1 || p == t2 || p == t3) return 1;
（原三个随机数）是可以AC的，但是上面那种写法才是理论正确的，所以多加了一个随机数吧tql
*/
bool s1 = Miller_Rabin(61, p), s2 = Miller_Rabin(24251, p);
bool s3 = Miller_Rabin(t1, p), s4 = Miller_Rabin(t2, p);
bool s5 = Miller_Rabin(t3, p), s6 = Miller_Rabin(t4, p);
return s1 && s2 & s3 && s4 && s5 && s6;
}难道你以为这就完了？你会发现这题的小数据就跟毒瘤一样，为什么呢？因为我们\(rand\)的随机数有可能大于当前的数\(p\)，这对答案是没有贡献的，那么我们需要点优化：
ll t1 = rand() % (p - 3) + 2这句话把随机的范围缩小到了\([2,p-1)\)。然后小数据就可以过了
但是，这还不是结尾，我们知道，4个随机数把合数判成质数的概率是\(\frac{1}{4}\)，我们需要缩小这个概率，那么就需要增加实验次数，于是我们可以把\(check\)和\(Miller \ Rabin\)写在一起，并增加一个“次数”\(repeat\)，在一定的时间复杂度内，我们进行多次实验。
下面是优化后的核心部分(\(repaet\)开到\(23\)貌似可以过了，当然，为了保险，可以开到\(50\))
inline bool Miller_Rabin(ll p, int repeat) {
if(p == 2 || p == 3) return true;
if(p % 2 == 0 || p == 1) return false;
//特判
register ll k = p - 1;
register int cnt = 0;
while(!(k & 1)) ++cnt, k >>= 1;
srand((unsigned)time(NULL));//种子
for(register int i = 0; i < repeat; i++) {
register ll test = rand() % (p - 3) + 2;
register ll x = qpow(test, k, p), y = 0;
for(register int j = 0; j < cnt; j++) { //二次探测逆过程
y = qmul(x, x, p);
if(y == 1 && x != 1 && x != (p - 1))
return false;
x = y;
}
if(y != 1) return false; //费马小定理
}
return true;
}所有的改进思路感谢巨佬@雪中舞
2. 快速幂 | 快速乘(相信这个我不用说了)
typedef long long ll;
ll quick_mul(ll x, ll y, ll mod) {
return (x*y-(ll)((long double)x/mod*y)*mod+mod)%mod;
}
ll quick_pow(ll x, ll m, ll mod) {
ll ans = 1;
while(m) {
if(m & 1)
ans = ans * x % mod;
x = x * x % mod;
m >>= 1;
}
return ans;
}3. \(gcd\)(这个我也不用说了吧)
两种办法：更相减损术 | 辗转相除法
下面给出辗转相除法的代码：
ll gcd(ll x,ll y) {
ll temp;
while(y) {
temp = x % y;
x = y;
y = temp;
}
return x;
}当然，这是最简单的\(gcd\)算法，我们还可以进行二进制优化：
详见约数相关
下面为Pollard-Rho 算法实现：
这里还要引入一个“生日悖论”，利用这个玄学的悖论，我们就可以快速质因数分解了。然后我们可以利用\(rand()\)函数进行随机，在一直找到一个不符合的数进入死循环时利用\(Floyd\)的判断环的方式，一个“正常倍速”，一个“二倍速”，两者“第一次相遇时”，说明“走完了一圈”，即进入死循环。
//参考代码
void Pollard_Rho(int x)
{
if(test(x)) {//素数测试
ans=max(x, ans);
return;
}
ll t1 = rand() % (x-1) + 1;
ll t2 = (qmul(t2, t2, x)+b) % x;
ll b = rand() % (x-1) + 1;
for(;t1 != t2;) {//Flord判断是否进入死循环
ll t = gcd(abs(t1-t2), x);
if(t != 1 && t != x)
{
Pollard_Rho(t),Pollard_Rho(x/t);
}
t1 = (qmul(t1, t1, x)+b) % x;
t2 = (qmul(t2, t2, x)+b) % x;
t2 = (qmul(t2, t2, x)+b) % x;//“两倍速”
}
}上述代码时间复杂度还没达到理想时间复杂度\(O(\sqrt [4]N)\)，主要瓶颈在于频繁的求\(gcd\)。
那么这里有个优化的玄学常数\(127\)，每取\(127\)样本才进行\(1\)次\(gcd\)，看起来是挺划算的？
下面是优化的核心部分代码(这真的是我手打的)：
inline void Pollard_Rho(ll x) {
if(check(x)) { //检查质数
ans = max(x, ans);
return;
}
ll s1 = rand() % (x - 1) + 1;
ll m = rand() % (x - 1) + 1, p = 1;
ll s2 = (qmul(s1, s1, x) + m) % x;
for(register ll i = 1; s1 != s2; i++) {
//边界：不进入死循环 即S1 != S2（Floyd判断环）
p = qmul(p, abs(s1 - s2), x);
if(p == 0) {
ll K = gcd(abs(s1 - s2), x);
if(K != 1 && K != x)
Pollard_Rho(K), Pollard_Rho(x / K);
return;
}
if(i % 127 == 0) {
p = gcd(p, x);
if(p != 1 && p != x) {
Pollard_Rho(p), Pollard_Rho(x / p);
return;
}
p = 1;
}
s1 = (qmul(s1, s1, x) + m) % x;
s2 = (qmul(s2, s2, x) + m) % x;
s2 = (qmul(s2, s2, x) + m) % x; //跑两倍
}
//以下部分：最后有可能一个环长不到127
p = gcd(p, x);
if(p != 1 && p != x) {
Pollard_Rho(p), Pollard_Rho(x / p);
return;
}
}对此有一定理解的读者可以尝试完成模板题：
P4718 【模板】Pollard-Rho算法 ，这也是上一章质数筛法的必做题目，下面给出参考代码：
\(Pollard-Rho\) 算法 \(Code:\)
#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<cmath>
#include<ctime>
using namespace std;
typedef long long ll;
ll ans,n;
inline ll qmul(ll x, ll y, ll mod) {
return (x*y-(ll)((long double)x/mod*y)*mod+mod)%mod;
}
inline ll qpow(ll x, ll m, ll mod) {
register ll ans = 1;
while(m) {
if(m & 1)
ans = qmul(ans, x, mod) % mod;
x = qmul(x, x, mod) % mod;
m >>= 1;
}
return ans;
}
inline bool FMLT(ll x, ll p) {
return qpow(x, p - 1, p) == 1;
}
inline bool Miller_Rabin(ll p, int repeat) {
if(p == 2 || p == 3) return true;
if(p % 2 == 0 || p == 1) return false;
register ll k = p - 1;
register int cnt = 0;
while(!(k & 1)) ++cnt, k >>= 1;
srand((unsigned)time(NULL));
for(register int i = 0; i < repeat; i++) {
register ll test = rand() % (p - 3) + 2;
register ll x = qpow(test, k, p), y = 0;
for(register int j = 0; j < cnt; j++) {
y = qmul(x, x, p);
if(y == 1 && x != 1 && x != (p - 1))
return false;
x = y;
}
if(y != 1) return false;
}
return true;
}
inline ll gcd(ll x, ll y) {
register ll temp;
while(y) {
temp = x % y;
x = y;
y = temp;
}
return x;
}
inline void Pollard_Rho(ll x) {
if(Miller_Rabin(x, 50)) {
ans = max(x, ans);
return;
}
register ll s1 = rand() % (x - 1) + 1;
register ll m = rand() % (x - 1) + 1, p = 1;
register ll s2 = (qmul(s1, s1, x) + m) % x;
for(register ll i = 1; s1 != s2; i++) {
p = qmul(p, abs(s1 - s2), x);
if(p == 0) {
register ll K = gcd(abs(s1 - s2), x);
if(K != 1 && K != x)
Pollard_Rho(K), Pollard_Rho(x / K);
return;
}
if(i % 127 == 0) {
p = gcd(p, x);
if(p != 1 && p != x) {
Pollard_Rho(p), Pollard_Rho(x / p);
return;
}
p = 1;
}
s1 = (qmul(s1, s1, x) + m) % x;
s2 = (qmul(s2, s2, x) + m) % x;
s2 = (qmul(s2, s2, x) + m) % x;
}
p = gcd(p, x);
if(p != 1 && p != x) {
Pollard_Rho(p), Pollard_Rho(x / p);
return;
}
}
namespace IN_OUT {
template <class T>
inline void read(T &x) {
static char c;
static bool op;
while(!isdigit(c = getchar()) && c != '-');
x = (op = c == '-')? 0: c-'0';
while(isdigit(c = getchar()))
x = x * 10 + c - '0';
if(op) x = ~x + 1;
}
template <class T>
inline void output(T x) {
register char buf[30], *tail = buf;
if(x == 0) putchar('0');
else {
if(x < 0) putchar('-'), x = ~x + 1;
for(; x; x /= 10) *++tail = x % 10 + '0';
for(; tail != buf; --tail)
putchar(*tail);
}
putchar('\n');
}
}
using namespace IN_OUT;
int main() {
read(n);
for(register int i = 1; i <= n; i++) {
register ll x;
read(x);
if(Miller_Rabin(x, 50))
printf("Prime\n");
else {
ans=0;
while(ans == 0)
Pollard_Rho(x);
output(ans);
}
}
return 0;
}下面是模题的评测记录，加了快读快输和一堆\(register\)，$ repeat$的值为50。
※章末练习
\(END\)
\(PS:\)
部分材料来自于李煜东《算法竞赛进阶指南》
\(Pollard-Pho\)算法 参考于:
质因数分解的题目较少，欢迎补充
标签：return,因数分解,ll,register,Pollard,Rho,浅谈,mod 来源： https://www.cnblogs.com/Ning-H/p/11567223.html
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。