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

Question

Question:

Given an array arr[] with N integers.

What is the maximum number of items that can be chosen from the array so that their GCD is greater than 1?

Example:

4
30 42 105 1

Answer: 3

Constransts

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

My take:

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);
    }
}

After many tries, my code still clearly got TLE, is there any more optimized solution for this problem?

Answer 1

Interesting task you have. I implemented two variants of solutions.

All algorithms that are used in my code are: Greatest Common Divisor (through Euclidean Algorithm), Binary Modular Exponentiation, Pollard Rho, Trial Division, Fermat Primality Test.

First variant called SolveCommon() iteratively finds all possible unique factors of all numbers by computing pairwise Greatest Common Divisor.

When all possible unique factors are found one can compute count of each unique factor inside each number. Finally maximal count for any factor will be final answer.

Second variant called SolveFactorize() finds all factor by doing factorization of each number using three algorithms: Pollard Rho, Trial Division, Fermat Primality Test.

Pollard-Rho factorization algorithm is quite fast, it has time complexity O(N^(1/4)), so for 64-bit number it will take around 2^16 iterations. To compare, Trial Division algorithm has complexity of O(N^(1/2)) which is square times slower than Pollard Rho. So in code below Pollard Rho can handle 64 bit inputs, although not very fast.

First variant SolveCommon() is much faster than second SolveFactorize(), especially if numbers are quite large, timings are provided in console output after following code.

Code below as an example provides test of random 100 numbers each 20 bit. 64 bit 1000 numbers are too large to handle by SolveFactorize() method, but SolveCommon() method solves 1000 64-bit numbers within 1-2 seconds.

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();    
}

Output:

common    325 time 0.061
factorize 325 time 0.005

Answer 2

I think you need to search among all possible prime numbers to find out which prime number can divide most element in the array.

Code:

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;
}

Note that in this way you cannot find out the value of the GCD, the code is based on that we dont need to know it.