[Pollard Rho模板]超快的大数质因数分解

发布于 2022-09-12  465 次阅读


使用试除法得到质因数的复杂度是\(O(\sqrt{n})\),如果数字较大比如\(10^{18}\)就束手无策了。

在上一篇文章有写到Miller Rabin可以在\(O(klog^{2}{n})\)的复杂度下判断一个数字是不是素数,这次要介绍的Pollard Rho(普拉德柔)算法就是结合Miller Rabin的快速素性检验来进行的一个复杂度大约为\(O(n^{\frac{1}{4}})\)来分解一个数字的所有质因数的算法。为什么说是大约呢,因为真不好算哈哈哈。

介绍 / Introduction

同样不深究原理,会用好模板就行了。

如果想要了解原理的请看:Pollard_Rho - ZenyZ - 博客园 (cnblogs.com)

简单来说就是随机找到一个因子(不一定是质因子),利用生日悖论来提高找到因子的概率,用倍增来提高找因子的速度,用累乘gcd来优化常数,不断地分治进行,如果到了一个质数就停止递归并将其加入到因子列表中。

例题:P4718 【模板】Pollard-Rho算法 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)

模板 / Templete

#include <bits/stdc++.h>
#define int long long
using namespace std;

//四个工具函数 
int qmi(int a, int b, int p) // 快速幂
{
    int res = 1;
    while (b)// 注意!中间结果可能溢出,需要使用__int128过渡
    {
        if (b & 1)res = (__int128)res * a % p; 
        a = (__int128)a * a % p, b >>= 1;
    }
    return res;
}
int getabs(int x){return x < 0 ? -x : x;}//绝对值 
int f(int t,int c,int x){return ((__int128)t * t + c) % x;}//随机化函数 
int gcd(int a,int b){return b == 0 ? a : gcd(b, a % b);}//求最大公因数

//判断x是否是质数 
bool isprime(int x)//Miller_Rabin
{
	//特判 
    if (x < 3)return x == 2;
    if (x % 2 == 0)return false;
    //用下面这几个数字判断不会出错 
    int A[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022}, d = x - 1, r = 0;
    while (d % 2 == 0)d >>= 1, ++r;
    for (auto a : A)
    {
        int v = qmi(a, d, x);
        if (v <= 1 || v == x - 1)continue;
        for (int i = 0; i < r; ++i)
        {
            v = (__int128)v * v % x;
            if (v == x - 1 && i != r - 1)
            {
                v = 1;
                break;
            }
           if (v == 1)return false;
        }
        if (v != 1)return false;
    }
    return true;
}

//找出x的一个因子 
int Pollard_Rho(int x)
{
	int s = 0, t = 0, c = 1LL * rand() % (x - 1) + 1;//[1, x - 1]
	for(int i = 1; true; i <<= 1, s = t)
	{
		int mul_sum = 1, g = 1;
		for(int j = 1;j <= i; ++ j)
		{
			
			t = f(t, c, x);
			mul_sum = (__int128)mul_sum * getabs(t - s) % x;
			if(j % 127 == 0)
			{
				g = gcd(mul_sum, x);
				if(g > 1)return g;
			}
		}
		g = gcd(mul_sum, x);
		if(g > 1)return g;
	}
	return -1;
}

//分治主体 
void find_factor(int x, vector<int> &v)
{
	if(isprime(x))return v.push_back(x), void();
	
	int yz = x;
	//一定要确保yz是小于x的,否则重新找 
	while(yz >= x)yz = Pollard_Rho(x); 
	find_factor(yz, v), find_factor(x / yz, v);
}

signed main()
{
	int N;cin >> N;
	for(int i = 1;i <= N; ++ i)
	{
		int x;cin >> x;
		if(isprime(x))cout << "Prime" << '\n';
		else
		{
			vector<int> v;
			find_factor(x, v);
			//找完之后所有的质数都在v这个数组中 
			//v的大小不会超过70,毕竟1e18也才2^63
			sort(v.begin(), v.end());
			cout << v.back() << '\n';
		}
		
	}
	return 0;
}