选择数组中的最大数量，使其 GCD > 1

Question

问题：

给定一个包含N个整数的数组arr[]。

可以从数组中选择的最大项目数是多少，以便其GCD大于1？

示例：

4
30 42 105 1
Answer: 3

Constransts

N <= 10^3
arr[i] <= 10^18

我的看法：

void solve(int i, int gcd, int chosen){
if(i > n){
maximize(res, chosen);
return;
}
solve(i+1, gcd, chosen);
if(gcd == -1) solve(i+1, arr[i], chosen+1);
else{
int newGcd = __gcd(gcd, arr[i]);
if(newGcd > 1) solve(i+1, newGcd, chosen+1);
}
}

经过多次尝试，我的代码仍然很明显地得到了TLE，有没有针对这个问题的更优化的解决方案？

Answer 1

您有一项有趣的任务。我实现了两种不同的解决方案。

在我的代码中使用的所有算法都是：最大公约数(通过欧几里得算法(、二进制模指数、Pollard-Rho、Trial Division、Fermat Primality Test。

第一个称为SolveCommon()的变体通过计算成对的最大公约数迭代地找到所有数字的所有可能的唯一因子。

当找到所有可能的唯一因子时，可以计算每个数字中每个唯一因子的计数。最后，任何因素的最大计数将是最终答案。

第二种变体称为SolveFactorize()，通过使用三种算法对每个数字进行因子分解来找到所有因子：Pollard-Rho、Trial Division和Fermat Primality Test。

Pollard-Rho因子分解算法速度很快，时间复杂度为O(N^(1/4))，因此对于64位数字，它需要大约2^16次迭代。相比之下，Trial Division算法的复杂度为O(N^(1/2))，比Pollard-Rho慢平方倍。因此，在PollardRow下面的代码中，Rho可以处理64位输入，尽管速度不是很快。

第一个变体SolveCommon()比第二个SolveFactorize()快得多，尤其是当数字很大时，在下面的代码之后，控制台输出中会提供定时。

下面的代码作为一个例子提供了随机100个数字的测试，每个数字20位。64位1000数字太大，SolveFactorize()方法无法处理，但SolveCommon()方法在1-2秒内解决1000个64位数字。

在线试用！

#include <cstdint>
#include <random>
#include <tuple>
#include <unordered_map>
#include <algorithm>
#include <set>
#include <iostream>
#include <chrono>
#include <cmath>
#include <map>
#define LN { std::cout << "LN " << __LINE__ << std::endl; }
using u64 = uint64_t;
using u128 = unsigned __int128;
static std::mt19937_64 rng{123}; //{std::random_device{}()};
auto CurTime() {
return std::chrono::high_resolution_clock::now();
}
static auto const gtb = CurTime();
double Time() {
return std::llround(std::chrono::duration_cast<
std::chrono::duration<double>>(CurTime() - gtb).count() * 1000) / 1000.0;
}
u64 PowMod(u64 a, u64 b, u64 const c) {
u64 r = 1;
while (b != 0) {
if (b & 1)
r = (u128(r) * a) % c;
a = (u128(a) * a) % c;
b >>= 1;
}
return r;
}
bool IsFermatPrp(u64 N, size_t ntrials = 24) {
// https://en.wikipedia.org/wiki/Fermat_primality_test
if (N <= 10)
return N == 2 || N == 3 || N == 5 || N == 7;
for (size_t trial = 0; trial < ntrials; ++trial)
if (PowMod(rng() % (N - 3) + 2, N - 1, N) != 1)
return false;
return true;
}
bool FactorTrialDivision(u64 N, std::vector<u64> & factors, u64 limit = u64(-1)) {
// https://en.wikipedia.org/wiki/Trial_division
if (N <= 1)
return true;
while ((N & 1) == 0) {
factors.push_back(2);
N >>= 1;
}
for (u64 d = 3; d <= limit && d * d <= N; d += 2)
while (N % d == 0) {
factors.push_back(d);
N /= d;
}
if (N > 1)
factors.push_back(N);
return N == 1;
}
u64 GCD(u64 a, u64 b) {
// https://en.wikipedia.org/wiki/Euclidean_algorithm
while (b != 0)
std::tie(a, b) = std::make_tuple(b, a % b);
return a;
}
bool FactorPollardRho(u64 N, std::vector<u64> & factors) {
// https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
auto f = [N](auto x) -> u64 { return (u128(x + 1) * (x + 1)) % N; };
auto DiffAbs = [](auto x, auto y){ return x >= y ? x - y : y - x; };
if (N <= 1)
return true;
if (IsFermatPrp(N)) {
factors.push_back(N);
return true;
}
for (size_t trial = 0; trial < 8; ++trial) {
u64 x = rng() % (N - 2) + 1;
size_t total_steps = 0;
for (size_t cycle = 1;; ++cycle) {
bool good = true;
u64 y = x;
for (u64 i = 0; i < (u64(1) << cycle); ++i) {
x = f(x);
++total_steps;
u64 const d = GCD(DiffAbs(x, y), N);
if (d > 1) {
if (d == N) {
good = false;
break;
}
//std::cout << N << ": " << d << ", " << total_steps << std::endl;
if (!FactorPollardRho(d, factors))
return false;
if (!FactorPollardRho(N / d, factors))
return false;
return true;
}
}
if (!good)
break;
}
}
factors.push_back(N);
return false;
}
void Factor(u64 N, std::vector<u64> & factors) {
if (N <= 1)
return;
if (1) {
FactorTrialDivision(N, factors, 1 << 8);
N = factors.back();
factors.pop_back();
}
FactorPollardRho(N, factors);
}
size_t SolveFactorize(std::vector<u64> const & nums) {
std::unordered_map<u64, size_t> cnts;
std::vector<u64> factors;
std::set<u64> unique_factors;
for (auto num: nums) {
factors.clear();
Factor(num, factors);
//std::cout << num << ": "; for (auto f: factors) std::cout << f << " "; std::cout << std::endl;
unique_factors.clear();
unique_factors.insert(factors.begin(), factors.end());
for (auto f: unique_factors)
++cnts[f];
}
size_t max_cnt = 0;
for (auto [_, cnt]: cnts)
max_cnt = std::max(max_cnt, cnt);
return max_cnt;
}
size_t SolveCommon(std::vector<u64> const & nums) {
size_t const K = nums.size();
std::set<u64> cmn(nums.begin(), nums.end()), cmn2, tcmn;
std::map<u64, bool> used;
cmn.erase(1);
while (true) {
cmn2.clear();
used.clear();
for (auto i = cmn.rbegin(); i != cmn.rend(); ++i) {
auto j = i;
++j;
for (; j != cmn.rend(); ++j) {
auto gcd = GCD(*i, *j);
if (gcd != 1) {
used[*i] = true;
used[*j] = true;
cmn2.insert(gcd);
cmn2.insert(*i / gcd);
cmn2.insert(*j / gcd);
break;
}
}
if (!used[*i])
tcmn.insert(*i);
}
cmn2.erase(1);
if (cmn2.empty())
break;
cmn = cmn2;
}
//for (auto c: cmn) std::cout << c << " "; std::cout << std::endl;
std::unordered_map<u64, size_t> cnts;
for (auto num: nums)
for (auto c: tcmn)
if (num % c == 0)
++cnts[c];
size_t max_cnt = 0;
for (auto [_, cnt]: cnts)
max_cnt = std::max(max_cnt, cnt);
return max_cnt;
}
void TestRandom() {
size_t const cnt_nums = 1000;

std::vector<u64> nums;
for (size_t i = 0; i < cnt_nums; ++i) {
nums.push_back((rng() & ((u64(1) << 20) - 1)) | 1);
//std::cout << nums.back() << " ";
}
//std::cout << std::endl;
{
auto tb = Time();
std::cout << "common    " << SolveCommon(nums) << " time " << (Time() - tb) << std::endl;
}
{
auto tb = Time();
std::cout << "factorize " << SolveFactorize(nums) << " time " << (Time() - tb) << std::endl;
}
}
int main() {
TestRandom();    
}

输出：

common    325 time 0.061
factorize 325 time 0.005

Answer 2

我认为您需要在所有可能的素数中搜索，以找出哪个素数可以划分数组中的大多数元素。

代码：

std::vector<int> primeLessEqualThanN(int N) {
std::vector<int> primes;    
for (int x = 2; x <= N; ++x) {
bool isPrime = true;
for (auto& p : primes) {
if (x % p == 0) {
isPrime = false;
break;
}
}
if (isPrime) primes.push_back(x);
}
return primes;
}
int maxNumberGCDGreaterThan1(int N, std::vector<int>& A) {
int A_MAX = *std::max_element(A.begin(), A.end());  // largest number in A
std::vector<int> primes = primeLessEqualThanN(std::sqrt(A_MAX)); 
int max_count = 0;
for (auto& p : primes) {
int count = 0;
for (auto& n : A)
if (n % p == 0)
count++;
max_count = count > max_count ? count : max_count;
}
return max_count;
}

请注意，通过这种方式，您无法找到GCD的值，代码基于此，我们不需要知道它。