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

Question

问题：

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

最多可以从数组中选择多少项，使其 GCD 大于 1？

示例：

4
30 42 105 1

Answer: 3

转换

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

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

代码：

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的值的，代码是基于我们不需要知道的。

Answer 2

你的任务很有趣。我实施了两种解决方案。

我的代码中使用的所有算法是：Greatest Common Divisor (through Euclidean Algorithm), Binary Modular Exponentiation, 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 慢了平方倍。所以在下面的代码中，Pollard Rho 可以处理 64 位输入，虽然不是很快。

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

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

Try it online!

#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

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

Choose maximum number in array so that their GCD is > 1

c++

algorithm

dynamic-programming