Первый шаг - использовать Сито Эратосфена для вычисления простых чисел до sqrt(10^9)
.Затем это сито можно использовать для быстрого поиска всех простых факторов любого числа, меньшего чем 10^9
(см. Функцию getPrimeFactors(...)
в примере кода ниже).
Далее, для каждого A[i]
с простыми коэффициентамиp0, p1, ..., pk
, мы вычисляем все возможные субпродукты X
- p0, p1, p0p1, p2, p0p2, p1p2, p0p1p2, p3, p0p3, ..., p0p1p2...pk
и считаем их на карте cntp[X]
.Фактически, карта cntp[X]
сообщает нам количество элементов A[i]
, делимое на X
, где X
- произведение простых чисел на степень 0 или 1. Так, например, для числа A[i] = 12
, главные факторы 2, 3
.Мы будем считать cntp[2]++
, cntp[3]++
и cntp[6]++
.
Наконец, для каждого B[j]
с простыми коэффициентами p0, p1, ..., pk
мы снова вычисляем все возможные субпродукты X
и используем Принцип включения-исключения для подсчета всех пар, не являющихся взаимно простыми C_j
(т. Е. Число A[i]
с, которые имеют хотя бы один простой множитель с B[j]
).Числа C_j
затем вычитаются из общего числа пар - N*N
, чтобы получить окончательный ответ.
Примечание: принцип включения-исключения выглядит следующим образом:
C_j = (cntp[p0] + cntp[p1] + ... + cntp[pk]) -
(cntp[p0p1] + cntp[p0p2] + ... + cntp[pk-1pk]) +
(cntp[p0p1p2] + cntp[p0p1p3] + ... + cntp[pk-2pk-1pk]) -
...
и объясняет тот факт, что в cntp[X]
и cntp[Y]
мы могли бы сосчитать одно и то же число A[i]
дважды, учитывая, что оно делится на X
и Y
.
возможная реализация алгоритма на C ++, которая дает те же результаты, что и простой алгоритм O (n ^ 2) с помощью OP:
// get prime factors of a using pre-generated sieve
std::vector<int> getPrimeFactors(int a, const std::vector<int> & primes) {
std::vector<int> f;
for (auto p : primes) {
if (p > a) break;
if (a % p == 0) {
f.push_back(p);
do {
a /= p;
} while (a % p == 0);
}
}
if (a > 1) f.push_back(a);
return f;
}
// find coprime pairs A_i and B_j
// A_i and B_i <= 1e9
void solution(const std::vector<int> & A, const std::vector<int> & B) {
// generate prime sieve
std::vector<int> primes;
primes.push_back(2);
for (int i = 3; i*i <= 1e9; ++i) {
bool isPrime = true;
for (auto p : primes) {
if (i % p == 0) {
isPrime = false;
break;
}
}
if (isPrime) {
primes.push_back(i);
}
}
int N = A.size();
struct Entry {
int n = 0;
int64_t p = 0;
};
// cntp[X] - number of times the product X can be expressed
// with prime factors of A_i
std::map<int64_t, int64_t> cntp;
for (int i = 0; i < N; i++) {
auto f = getPrimeFactors(A[i], primes);
// count possible products using non-repeating prime factors of A_i
std::vector<Entry> x;
x.push_back({ 0, 1 });
for (auto p : f) {
int k = x.size();
for (int i = 0; i < k; ++i) {
int nn = x[i].n + 1;
int64_t pp = x[i].p*p;
++cntp[pp];
x.push_back({ nn, pp });
}
}
}
// use Inclusion–exclusion principle to count non-coprime pairs
// and subtract them from the total number of prairs N*N
int64_t cnt = N; cnt *= N;
for (int i = 0; i < N; i++) {
auto f = getPrimeFactors(B[i], primes);
std::vector<Entry> x;
x.push_back({ 0, 1 });
for (auto p : f) {
int k = x.size();
for (int i = 0; i < k; ++i) {
int nn = x[i].n + 1;
int64_t pp = x[i].p*p;
x.push_back({ nn, pp });
if (nn % 2 == 1) {
cnt -= cntp[pp];
} else {
cnt += cntp[pp];
}
}
}
}
printf("cnt = %d\n", (int) cnt);
}
Живой пример
Не могуОцените сложность аналитически, но вот некоторые результаты профилирования на моем ноутбуке для различных N
и равномерно случайных A[i]
и B[j]
:
For N = 1e2, takes ~0.02 sec
For N = 1e3, takes ~0.05 sec
For N = 1e4, takes ~0.38 sec
For N = 1e5, takes ~3.80 sec
Для сравнения, подход O (n ^ 2)занимает:
For N = 1e2, takes ~0.00 sec
For N = 1e3, takes ~0.15 sec
For N = 1e4, takes ~15.1 sec
For N = 1e5, takes too long, didn't wait to finish