浅谈质因数分解
->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\)算法 参考于:
质因数分解的题目较少，欢迎补充