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

Time:12-15

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?

CodePudding user response:

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.

CodePudding user response:

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

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.

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 32 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.

If my SolveCommon() solution gives timeout on your largest test, then tell me in comments, I have few ideas how to make it faster.

Try it online!

#include <cstdint>
#include <random>
#include <tuple>
#include <unordered_map>
#include <algorithm>
#include <set>
#include <iostream>
#include <chrono>
#include <cmath>

#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 c) {
    u64 r = 1;
    while (c != 0) {
        if (c & 1)
            r = (u128(r) * a) % c;
        a = (u128(a) * a) % c;
        c >>= 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) {
    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's_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;
    cmn.erase(1);
    while (true) {
        cmn2.clear();
        for (auto i = cmn.rbegin(); i != cmn.rend();   i) {
            auto j = i;
              j;
            bool found = false;
            for (; j != cmn.rend();   j) {
                auto gcd = GCD(*i, *j);
                if (gcd != 1) {
                    found = true;
                    cmn2.insert(gcd);
                    cmn2.insert(*i / gcd);
                    cmn2.insert(*j / gcd);
                    break;
                }
            }
            if (!found)
                cmn2.insert(*i);
        }
        cmn2.erase(1);
        if (cmn.size() == cmn2.size() && cmn == cmn2)
            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: cmn)
            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 = 100;
    
    std::vector<u64> nums;
    for (size_t i = 0; i < cnt_nums;   i) {
        nums.push_back((rng() & ((u64(1) << 32) - 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    38 time 0.007
factorize 38 time 4.347

CodePudding user response:

For rather small limitation N<1000 you can calculate GCD for all item pairs - at most 5*10^5 calculations.

Now items form undirected graph where edge between nodes a,b exists when GCD(a,b)>1

Your task is to build (or just find size) of maximal clique (subset of nodes, all adjacent to each other).

Sadly, this is NP-Hard problem (requires exponential time), and using both brute-force and advanced algorithms (like Bron–Kerbosch one) looks unsolvable for N=1000.

  • Related