PROBLEM LINK:
Practice
Contest: Division 1
Contest: Division 2
Contest: Division 3
Contest: Division 4
Author: gudegude
Testers: IceKnight1093, tejas10p
Editorialist: IceKnight1093
DIFFICULTY:
2649
PREREQUISITES:
PROBLEM:
Given an array A, find the smallest positive integer d such that d doesn’t divide |A_i - A_j| for any 1 \leq i \lt j \leq N.
EXPLANATION:
Let M = \max(A) be the maximum element of the array.
Note that any difference between array elements will also be \leq M.
So, if we were able to compute every pair of differences, we could solve the problem fairly easily, almost brute-force:
- Iterate d from 1 to M+1
- Once d is fixed, check whether d, 2d, 3d, \ldots appear as a difference of array elements.
- If none of them appear the answer is d; otherwise increment d by 1 and try again.
Assuming we can check whether x appears as a difference in \mathcal{O}(1), the complexity of the above algorithm is \displaystyle \mathcal{O}\left(\frac{M}{1} + \frac{M}{2} + \frac{M}{3} + \ldots\right) = \mathcal{O}(M\log M) (it’s the harmonic series multiplied by M).
Finding all differences
The task of finding all pairwise differences among array elements (or between two different arrays) is in fact a rather standard application of FFT, the fast fourier transform.
The idea is to use polynomials, as follows:
Consider two polynomials p(x) = p_0 + p_1x + p_2x^2 + \ldots and q(x) = q_0 + q_1x + q_2x^2 +\ldots, where:
- p_i is the number of times i occurs in A.
- q_i is the number of times M-i occurs in A (recall that M = \max(A))
p and q are degree-(M+1) polynomials.
Let r(x) = p(x) \cdot q(x) be their product.
Note that for any x and y, such that x \leq y \leq M, exactly p_y \cdot q_{M-x} is added to r_{M+y-x}.
In particular, for d \geq 0, r_{M+d} is exactly the number of pairs whose difference is d, which is exactly what we want!
r is the product of two polynomials, which can be computed in \mathcal{O}(M\log M) using FFT.
TIME COMPLEXITY:
\mathcal{O}(N + M \log M) per testcase, where M = \max(A).
CODE:
Tester's code (C++)
#include <bits/stdc++.h>
#define ll long long int
using namespace std;
namespace FFT {
using float_t = double;
struct cd {
float_t x, y;
cd(float_t x = 0, float_t y = 0): x(x), y(y) {}
friend ostream& operator<<(ostream& os, const cd& c) {
return os << '(' << c.x << ", " << c.y << ')';
}
float_t real() const { return x; }
float_t imag() const { return y; }
cd& operator+=(const cd& o) { x += o.x, y += o.y; return *this; }
cd operator+(const cd& o) const { cd ret = *this; ret += o; return ret; }
cd& operator-=(const cd& o) { x -= o.x, y -= o.y; return *this; }
cd operator-(const cd& o) const { cd ret = *this; ret -= o; return ret; }
cd& operator*=(const cd& o) { tie(x, y) = make_pair(x * o.x - y * o.y, x * o.y + y * o.x); return *this; }
cd operator*(const cd& o) const { cd ret = *this; ret *= o; return ret; }
cd& operator/=(const float_t& r) { x /= r, y /= r; return *this; }
cd operator/(const float_t& r) const { cd ret = *this; ret /= r; return ret; }
};
const float_t PI = acos(float_t(-1));
const int FFT_CUTOFF = 150;
vector<cd> roots = {{0, 0}, {1, 0}}, inv_roots = {{0, 0}, {1, 0}};
vector<int> bit_order;
void prepare_roots(int n) {
if(roots.size() >= n) return;
int len = roots.size();
roots.resize(n); inv_roots.resize(n);
while(n > len) {
const cd w(cos(PI / len), sin(PI / len)), inv_w(w.real(), -w.imag());
for(int i = len >> 1; i < len; i++) {
roots[i << 1] = roots[i];
roots[i << 1 | 1] = roots[i] * w;
inv_roots[i << 1] = inv_roots[i];
inv_roots[i << 1 | 1] = inv_roots[i] * inv_w;
} len <<= 1;
}
}
void reorder_bits(int n, vector<cd>& a) {
if(bit_order.size() != n) {
bit_order.assign(n, 0);
int len = __builtin_ctz(n);
for(int i = 0; i < n; i++)
bit_order[i] = (bit_order[i >> 1] >> 1) + ((i & 1) << len-1);
}
for(int i = 0; i < n; i++)
if(i > bit_order[i]) swap(a[i], a[bit_order[i]]);
}
void fft(int n, vector<cd>& a, bool invert = false) {
prepare_roots(n); reorder_bits(n, a);
const auto& ws = invert? inv_roots : roots;
for(int len = 1; len < n; len <<= 1)
for(int i = 0; i < n; i += len << 1)
for(int j = 0; j < len; j++) {
const cd& even = a[i + j];
cd odd = a[i + len + j] * ws[len + j];
a[i + len + j] = even - odd; a[i + j] = even + odd;
}
if(invert)
for(auto& x: a) x /= n;
}
template<typename T>
vector<T> multiply(const vector<T>& a, const vector<T>& b, bool trim = false) {
int n = a.size(), m = b.size();
if(n == 0 or m == 0)
return vector<T>{0};
if(min(n, m) < FFT_CUTOFF) {
vector<T> res(n + m - 1);
for(int i = 0; i < n; i++)
for(int j = 0; j < m; j++)
res[i + j] += a[i] * b[j];
if(trim) {
while(!res.empty() && res.back() == 0)
res.pop_back();
}
return res;
}
int N = [](int x) { while(x & x-1) x = (x | x-1) + 1; return x; }(n + m - 1);
vector<cd> fa(a.begin(), a.end()), fb(b.begin(), b.end());
fa.resize(N); fb.resize(N);
bool equal = n == m and a == b;
fft(N, fa);
if(equal) fb = fa;
else fft(N, fb);
for(int i = 0; i < N; i++)
fa[i] *= fb[i];
fft(N, fa, true);
fa.resize(n + m - 1);
vector<T> res(n + m - 1);
for(int i = 0; i < n + m - 1; i++)
res[i] = is_integral<T>::value? round(fa[i].real()) : fa[i].real();
if(trim) {
while(!res.empty() && res.back() == 0)
res.pop_back();
}
return res;
}
template<typename T>
T expo(T A, int64_t B) {
T res{1}; while(B) {
if(B & 1) res = res * A;
B >>= 1; A = A * A;
} return res;
}
template<typename T>
vector<T> expo(const vector<T>& a, int e, bool trim = false) {
int n = a.size();
int N = [](int x) { while(x & x-1) x = (x | x-1) + 1; return x; }((n-1) * e + 1);
vector<cd> fa(a.begin(), a.end()); fa.resize(N);
fft(N, fa);
for(int i = 0; i < N; i++)
fa[i] = expo(fa[i], e);
fft(N, fa, true);
fa.resize((n-1) * e + 1);
vector<T> res((n-1) * e + 1);
for(int i = 0; i < res.size(); i++)
res[i] = is_integral<T>::value? round(fa[i].real()) : fa[i].real();
if(trim) {
while(!res.empty() && res.back() == 0)
res.pop_back();
}
return res;
}
} // namespace FFT
int main() {
int t;
cin >> t;
while(t--) {
int n;
cin >> n;
vector<int> v(n);
for(int i = 0; i < n; i++) cin >> v[i];
sort(v.begin(), v.end());
ll NMAX = 2*v[n - 1] + 1;
vector<ll> a(NMAX, 0), b(NMAX, 0);
for(int i = 0; i < n; i++) a[v[n - 1] + v[i]] = 1, b[v[n - 1] - v[i]] = 1;
vector<ll> c = FFT::multiply<ll>(a, b);
for(int i = 1; ; i++) {
bool good = 1;
for(int j = i; j <= 2*v[n - 1]; j+=i) {
if(c[2*v[n - 1] + j] > 0.5) good = 0;
}
if(good) {
cout << i << "\n";
break;
}
}
}
return 0;
}
Editorialist's code (C++)
#include "bits/stdc++.h"
// #pragma GCC optimize("O3,unroll-loops")
// #pragma GCC target("avx2,bmi,bmi2,lzcnt,popcnt")
using namespace std;
using ll = long long int;
mt19937_64 rng(chrono::high_resolution_clock::now().time_since_epoch().count());
const long double PI = 2*acos(0.0);
void fft(vector<complex<double>>& a) {
int n = size(a), L = 31 - __builtin_clz(n);
static vector<complex<long double>> R(2, 1);
static vector<complex<double>> rt(2, 1); // (^ 10% faster if double)
for (static int k = 2; k < n; k *= 2) {
R.resize(n); rt.resize(n);
const long double PI = 2*acos(0.0);
auto x = polar(1.0L, PI / k); // M_PI, lower-case L
for (int i = k; i < 2*k; ++i) rt[i] = R[i] = i&1 ? R[i/2] * x : R[i/2];
}
vector<int> rev(n);
for (int i = 0; i < n; ++i) rev[i] = (rev[i / 2] | (i & 1) << L) / 2;
for (int i = 0; i < n; ++i) if (i < rev[i]) swap(a[i], a[rev[i]]);
for (int k = 1; k < n; k *= 2)
for (int i = 0; i < n; i += 2 * k) for (int j = 0; j < k; ++j) {
complex<double> z = rt[j+k] * a[i+j+k]; // (25% faster if hand-rolled) /// include-line
a[i + j + k] = a[i + j] - z;
a[i + j] += z;
}
}
vector<double> conv(const vector<double>& a, const vector<double>& b) {
if (a.empty() || b.empty()) return {};
vector<double> res(size(a) + size(b) - 1);
int L = 32 - __builtin_clz(size(res)), n = 1 << L;
vector<complex<double>> in(n), out(n);
copy(begin(a), end(a), begin(in));
for (int i = 0; i < (int)size(b); ++i) in[i].imag(b[i]);
fft(in);
for (auto &x : in) x *= x;
for (int i = 0; i < n; ++i) out[i] = in[-i & (n - 1)] - conj(in[i]);
fft(out);
for (int i = 0; i < (int)size(res); ++i) res[i] = imag(out[i]) / (4 * n);
return res;
}
int main()
{
ios::sync_with_stdio(false); cin.tie(0);
int t; cin >> t;
while (t--) {
int n; cin >> n;
vector<int> a(n);
for (int &x : a) cin >> x;
int mx = *max_element(begin(a), end(a));
vector<double> forward(mx+1), backward(mx+1);
for (int x : a) {
++forward[x];
++backward[mx-x];
}
auto prod = conv(forward, backward);
// x and y in a -> add 1 to forward[x] * backward[mx - y] -> prod[mx + (x - y)]
int ans;
for (ans = 1; ; ++ans) {
bool found = false;
for (int i = ans; !found; i += ans) {
if (mx + i >= prod.size()) break;
found |= prod[mx + i] > 0.5;
}
if (!found) break;
}
cout << ans << '\n';
}
}