This solution works in time O(n + h log h), where h is the maximum number in the array. Let's solve a harder problem: for each x <= h, count the number d[x] of unordered pairs (i, j) such that 0 <= i, j < n and GCD(a[i], a[j]) = x. To solve your problem, just find the smallest x such that d[x] is positive. Also note that counting ordered pairs (i, j) with i = j does not influence the solution. The solution uses Mobius inversion - basically a variation of Inclusion-Exclusion for divisors of integers.
Mobius inversion can be used to solve the following problem: you need to find an array y, but you're given an array z such that z[k] = y[k] + y[2*k] + y[3*k] + .... Surprisingly, it works in-place and it's just three lines of code!
This is exactly what we need, first we'll find the number of ordered pairs (i, j) such that d[x] divides GCD(a[i], a[j]), but we need the number of ordered pairs (i, j) such that d[x] is GCD(a[i], a[j]).
#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;
typedef long long ll;
int main() {
int n, h = 0;
cin >> n;
vector<int> a(n);
for (int& x : a) {
cin >> x;
h = max(h, x);
}
h++;
vector<ll> c(h), d(h);
for (int x : a)
c[x]++;
for (int i=1; i<h; i++)
for (int j=i; j<h; j+=i)
d[i] += c[j];
// now, d[x] = no. of indices i such that x divides a[i]
for (int i=1; i<h; i++)
d[i] *= d[i];
// now, d[x] = number of pairs of indices (i, j) such that
// x divides a[i] and a[j] (just square the previous d)
// equiv. x divides GCD(a[i], a[j])
// apply Mobius inversion to get proper values of d
for (int i=h-1; i>0; i--)
for (int j=2*i; j<h; j+=i)
d[i] -= d[j];
for (int i=1; i<h; i++) {
if (d[i]) {
cout << i << '\n';
return 0;
}
}
}