质因数分解最大为 10^18 的最快方法
Fastest way to prime factorise a number up to 10^18
给定一个数字 1 <= n <= 10^18，我怎样才能以最少的时间复杂度分解它？
互联网上有很多帖子介绍如何找到主要因素，但其中 none（至少从我所见）陈述了它们的好处，比如在特定情况下。
除了 Eratosthenes 的筛法之外，我还使用 Pollard 的 rho 算法：
- 用筛法，找出前107个数中的所有质数，然后用这些质数尽可能多地除以
n
- 现在使用 Pollard 的 rho 算法尝试找到其余素数，直到 n 等于 1。
我的实现：
#include <iostream>
#include <vector>
#include <cstdio>
#include <ctime>
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <string>
using namespace std;
typedef unsigned long long ull;
typedef long double ld;
typedef pair <ull, int> pui;
#define x first
#define y second
#define mp make_pair
bool prime[10000005];
vector <ull> p;
void initprime(){
prime[2] = 1;
for(int i = 3 ; i < 10000005 ; i += 2){
prime[i] = 1;
}
for(int i = 3 ; i * i < 10000005 ; i += 2){
if(prime[i]){
for(int j = i * i ; j < 10000005 ; j += 2 * i){
prime[j] = 0;
}
}
}
for(int i = 0 ; i < 10000005 ; ++i){
if(prime[i]){
p.push_back((ull)i);
}
}
}
ull modularpow(ull base, ull exp, ull mod){
ull ret = 1;
while(exp){
if(exp & 1){
ret = (ret * base) % mod;
}
exp >>= 1;
base = (base * base) % mod;
}
return ret;
}
ull gcd(ull x, ull y){
while(y){
ull temp = y;
y = x % y;
x = temp;
}
return x;
}
ull pollardrho(ull n){
srand(time(NULL));
if(n == 1)
return n;
ull x = (rand() % (n - 2)) + 2;
ull y = x;
ull c = (rand() % (n - 1)) + 1;
ull d = 1;
while(d == 1){
x = (modularpow(x, 2, n) + c + n) % n;
y = (modularpow(y, 2, n) + c + n) % n;
y = (modularpow(y, 2, n) + c + n) % n;
d = gcd(abs(x - y), n);
if(d == n){
return pollardrho(n);
}
}
return d;
}
int main ()
{
ios_base::sync_with_stdio(false);
cin.tie(0);
initprime();
ull n;
cin >> n;
ull c = n;
vector <pui> o;
for(vector <ull>::iterator i = p.begin() ; i != p.end() ; ++i){
ull t = *i;
if(!(n % t)){
o.push_back(mp(t, 0));
}
while(!(n % t)){
n /= t;
o[o.size() - 1].y++;
}
}
while(n > 1){
ull u = pollardrho(n);
o.push_back(mp(u, 0));
while(!(n % u)){
n /= u;
o[o.size() - 1].y++;
}
if(n < 10000005){
if(prime[n]){
o.push_back(mp(n, 1));
}
}
}
return 0;
}
有没有更快的方法来分解这些数字？如果可能，请在源代码中说明原因。
现代处理器上 64 位输入的最快解决方案是进行少量试验除法（数量会有所不同，但通常低于 100），然后是 Pollard 的 Rho。您将需要使用 Miller-Rabin 或 BPSW 进行良好的确定性素数测试，以及处理多个因素的基础设施（例如，如果将复合材料拆分为更多复合材料）。对于 32 位，您可以进一步优化这些内容。
您需要一个快速的 mulmod，因为它是 Pollard 的 Rho、Miller-Rabin 和 Lucas 检验的核心。理想情况下，这是作为一个小的汇编程序片段完成的。
分解任何 64 位输入的时间应低于 1 毫秒。在 50 位以下明显更快。
如 Ben Buhrow 的 spBrent 实现所示，Brent 1980 年论文中的算法 P2'' 似乎与我所知道的其他实现一样快。它使用 Brent 改进的周期查找以及通过必要的添加回溯延迟 GCD 的有用技巧。
有关各种解决方案的一些杂乱细节和基准测试，请参阅 this thread on Mersenneforum。我有许多不同大小的这些和其他实现的基准，但还没有发布任何东西（部分原因是有很多方法可以查看数据）。
由此产生的真正有趣的事情之一是 SQUFOF，多年来被认为是 64 位范围高端的更好解决方案，不再具有竞争力。 SQUFOF 的优势在于只需要一个快速的完美方形检测器即可获得最佳速度，而不必在 asm 中才能非常快。
方法
假设您有一个数 n，最大为 1018，并且您想对它进行质因数分解。因为这个数字可以小到 1018，一直以来它既可以是素数也可以是合数，这就是我的方法 -
- 使用miller rabin primality testing，确保数字是合数。
- 使用不超过 106 的素数因式分解
n，可以使用 sieve of Eratosthenes. 进行计算
现在 n 的更新值是这样的，它只有大于 106 的质因数，因为 n 的值仍然可以是大到 1018，我们得出结论，这个数要么是素数，要么恰好有两个素数因子（不一定不同）。
- 运行 Miller Rabin 再次确保数字不是质数。
- 使用Pollard rho algorithm得到一个质因数。
你现在已经完成了因式分解。
让我们看看上述方法的时间复杂度：
- 米勒拉宾拿下
O(log n)
- Eratosthenes 筛法需要
O(n*log n)
- 我分享的 Pollard rho 的实现需要
O(n^0.25)
时间复杂度
第2步花费的最大时间等于O(10^7)，这又是上述算法的复杂度。这意味着您可以在一秒钟内找到几乎所有编程语言的因式分解。
Space 复杂度
Space仅在执行sieve的第2步中使用，等于O(10^6)。同样，非常实用。
实施
Complete Code 在 C++14 中实施。该代码有一个隐藏的错误。您可以在下一节中展示它，也可以跳到挑战 ;)
代码中的错误
在 line 105 中，迭代直到 i<=np。否则，您可能会错过 prime[np]=999983 是素数的情况
挑战
给我一个 n 的值，如果有的话，共享代码导致错误的质因数分解。
奖金
存在多少个 n 这样的值？
暗示
对于这样的 n 值，Line 119 中的断言可能会失败。
解决方案
让我们打电话给 P=999983。 n = p*q*r 形式的所有数字，其中 p、q、r 是素数 >= P，使得其中至少一个等于 P 将导致错误的素数分解。
奖金解决方案
有正好有四个这样的数字：{P03, P02P1, P0 2P2, P0P12}，其中 P0 = P = 999983，P1 = next_prime(P0 ) = 1000003, P2 = next_prime(P1) = 1000033.
给定一个数字 1 <= n <= 10^18，我怎样才能以最少的时间复杂度分解它？
互联网上有很多帖子介绍如何找到主要因素，但其中 none（至少从我所见）陈述了它们的好处，比如在特定情况下。
除了 Eratosthenes 的筛法之外，我还使用 Pollard 的 rho 算法：
- 用筛法，找出前107个数中的所有质数，然后用这些质数尽可能多地除以
n - 现在使用 Pollard 的 rho 算法尝试找到其余素数，直到 n 等于 1。
我的实现：
#include <iostream>
#include <vector>
#include <cstdio>
#include <ctime>
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <string>
using namespace std;
typedef unsigned long long ull;
typedef long double ld;
typedef pair <ull, int> pui;
#define x first
#define y second
#define mp make_pair
bool prime[10000005];
vector <ull> p;
void initprime(){
prime[2] = 1;
for(int i = 3 ; i < 10000005 ; i += 2){
prime[i] = 1;
}
for(int i = 3 ; i * i < 10000005 ; i += 2){
if(prime[i]){
for(int j = i * i ; j < 10000005 ; j += 2 * i){
prime[j] = 0;
}
}
}
for(int i = 0 ; i < 10000005 ; ++i){
if(prime[i]){
p.push_back((ull)i);
}
}
}
ull modularpow(ull base, ull exp, ull mod){
ull ret = 1;
while(exp){
if(exp & 1){
ret = (ret * base) % mod;
}
exp >>= 1;
base = (base * base) % mod;
}
return ret;
}
ull gcd(ull x, ull y){
while(y){
ull temp = y;
y = x % y;
x = temp;
}
return x;
}
ull pollardrho(ull n){
srand(time(NULL));
if(n == 1)
return n;
ull x = (rand() % (n - 2)) + 2;
ull y = x;
ull c = (rand() % (n - 1)) + 1;
ull d = 1;
while(d == 1){
x = (modularpow(x, 2, n) + c + n) % n;
y = (modularpow(y, 2, n) + c + n) % n;
y = (modularpow(y, 2, n) + c + n) % n;
d = gcd(abs(x - y), n);
if(d == n){
return pollardrho(n);
}
}
return d;
}
int main ()
{
ios_base::sync_with_stdio(false);
cin.tie(0);
initprime();
ull n;
cin >> n;
ull c = n;
vector <pui> o;
for(vector <ull>::iterator i = p.begin() ; i != p.end() ; ++i){
ull t = *i;
if(!(n % t)){
o.push_back(mp(t, 0));
}
while(!(n % t)){
n /= t;
o[o.size() - 1].y++;
}
}
while(n > 1){
ull u = pollardrho(n);
o.push_back(mp(u, 0));
while(!(n % u)){
n /= u;
o[o.size() - 1].y++;
}
if(n < 10000005){
if(prime[n]){
o.push_back(mp(n, 1));
}
}
}
return 0;
}
有没有更快的方法来分解这些数字？如果可能，请在源代码中说明原因。
现代处理器上 64 位输入的最快解决方案是进行少量试验除法（数量会有所不同，但通常低于 100），然后是 Pollard 的 Rho。您将需要使用 Miller-Rabin 或 BPSW 进行良好的确定性素数测试，以及处理多个因素的基础设施（例如，如果将复合材料拆分为更多复合材料）。对于 32 位，您可以进一步优化这些内容。
您需要一个快速的 mulmod，因为它是 Pollard 的 Rho、Miller-Rabin 和 Lucas 检验的核心。理想情况下，这是作为一个小的汇编程序片段完成的。
分解任何 64 位输入的时间应低于 1 毫秒。在 50 位以下明显更快。
如 Ben Buhrow 的 spBrent 实现所示，Brent 1980 年论文中的算法 P2'' 似乎与我所知道的其他实现一样快。它使用 Brent 改进的周期查找以及通过必要的添加回溯延迟 GCD 的有用技巧。
有关各种解决方案的一些杂乱细节和基准测试，请参阅 this thread on Mersenneforum。我有许多不同大小的这些和其他实现的基准，但还没有发布任何东西（部分原因是有很多方法可以查看数据）。
由此产生的真正有趣的事情之一是 SQUFOF，多年来被认为是 64 位范围高端的更好解决方案，不再具有竞争力。 SQUFOF 的优势在于只需要一个快速的完美方形检测器即可获得最佳速度，而不必在 asm 中才能非常快。
方法
假设您有一个数 n，最大为 1018，并且您想对它进行质因数分解。因为这个数字可以小到 1018，一直以来它既可以是素数也可以是合数，这就是我的方法 -
- 使用miller rabin primality testing，确保数字是合数。
- 使用不超过 106 的素数因式分解
n，可以使用 sieve of Eratosthenes. 进行计算
现在 n 的更新值是这样的，它只有大于 106 的质因数，因为 n 的值仍然可以是大到 1018，我们得出结论，这个数要么是素数，要么恰好有两个素数因子（不一定不同）。
- 运行 Miller Rabin 再次确保数字不是质数。
- 使用Pollard rho algorithm得到一个质因数。
你现在已经完成了因式分解。
让我们看看上述方法的时间复杂度：
- 米勒拉宾拿下
O(log n) - Eratosthenes 筛法需要
O(n*log n) - 我分享的 Pollard rho 的实现需要
O(n^0.25)
时间复杂度
第2步花费的最大时间等于O(10^7)，这又是上述算法的复杂度。这意味着您可以在一秒钟内找到几乎所有编程语言的因式分解。
Space 复杂度
Space仅在执行sieve的第2步中使用，等于O(10^6)。同样，非常实用。
实施
Complete Code 在 C++14 中实施。该代码有一个隐藏的错误。您可以在下一节中展示它，也可以跳到挑战 ;)
代码中的错误
在
line 105中，迭代直到i<=np。否则，您可能会错过prime[np]=999983是素数的情况
挑战
给我一个 n 的值，如果有的话，共享代码导致错误的质因数分解。
奖金
存在多少个 n 这样的值？
解决方案对于这样的 n 值，
Line 119中的断言可能会失败。
奖金解决方案让我们打电话给
P=999983。n = p*q*r形式的所有数字，其中 p、q、r 是素数 >=P，使得其中至少一个等于P将导致错误的素数分解。
有正好有四个这样的数字：{P03, P02P1, P0 2P2, P0P12}，其中 P0 = P = 999983，P1 = next_prime(P0 ) = 1000003, P2 = next_prime(P1) = 1000033.