//Konrad Paluszek,University of Warsaw(former XIV LO Staszic)
//#STAY AT HOME
#ifndef LOCAL
#pragma GCC optimize("O3")
#endif
#define TIME (chrono::steady_clock::now().time_since_epoch().count())
#include<bits/stdc++.h>
#include<ext/pb_ds/assoc_container.hpp>
#include<ext/pb_ds/tree_policy.hpp>
#define xfm(a,b)a##b
#define xwb(a,b)xfm(a,b)
#define _ xwb(nvj,__LINE__)
#define __ xwb(kjy,__LINE__)
#define ___ xwb(cjm,__LINE__)
#define REP(i,n)for(urs(n)i=0;i<(n);++i)
#define PER(r...)for(bool _=1;_||next_permutation(r);_=false)
#define FS(r)r.first,r.second
#define ALL(r)(r).begin(),(r).end()
#define M0(r) memset(r,0,sizeof(r))
#define sim template<class c
#define forbits(i,m)if(m)for(urs(m)i=ctz(m),i##nvj=m;i##nvj;i##nvj^=((urs(m))1<<i),i=i##nvj?ctz(i##nvj):0)
#define fordbits(i,m)if(m)for(urs(m)i=8*sizeof(m)-clz(m)-1,i##nxd=m;i##nxd;i##nxd^=((urs(m))1<<i),i=i##nxd?8*sizeof(m)-clz(i##nxd)-1:0)
#define ksets(t, m, k, n)for(t m=(((t)1<<(k))-1);m<((t)1<<(n));m=nux(m))
#define urs(r...)typename decay<decltype(r)>::type
#define hur(f,g)sim>int f(c a){if(sizeof(c)==8)return g##ll(a);return g(a);}
using namespace __gnu_pbds;using namespace std;using ll=long long;using ld=long double;using ull=unsigned long long;using vi=vector<int>;using vll=vector<ll>;using pii=pair<int,int>;using pll=pair<ll,ll>;using vpii=vector<pii>;using unt=unsigned int;sim>using min_queue=priority_queue<c,vector<c>,greater<>>;sim,class b,class cmp=less<c>>using ordered_map=tree<c,b,cmp,rb_tree_tag,tree_order_statistics_node_update>;sim,class cmp=less<c>>using ordered_set=ordered_map<c,null_type,cmp>;hur(popc,__builtin_popcount)hur(ctz,__builtin_ctz)hur(clz,__builtin_clz)sim,class N>bool mini(c&o,const N&h){if(o>h)return o=h,1;return 0;}sim,class N>bool maxi(c&o,const N&h){if(o<h)return o=h,1;return 0;}
#ifdef LOCAL
#include </home/konrad/cp/headers/debuglib.hpp>
#else
#define loc(...)
#define onl(r...)r
#define debug(...)
#define print_stack(...)
#define mark_stack(...)
#define set_pre(...)
#define reg_it(...)
#define def_op(...) struct _{};
#define mask_set(...)
#define exit my_exit
void __resources();
void my_exit(int x) {cout.flush();
#ifdef LOCAL2
__resources();
#endif
_Exit(x);}
#endif
#define next nexT
#define prev preV
#define tree trEE
#define left lefT
#define right righT
#define div diV
#define y1 y_1
#define pow don't'
ull mix(ull o){o+=0x9e3779b97f4a7c15;o=(o^(o>>30))*0xbf58476d1ce4e5b9;o=(o^(o>>27))*0x94d049bb133111eb;return o^(o>>31);}ull SALT=0x7a14a4b0881ebf9,tqu=0x7a14a4b0881ebf9;ull my_rand(){return tqu=mix(tqu);}void my_srand(ull x){SALT=tqu=x;}const int inf=1023400000;const ll llinf=1234567890000000000ll;ll fix(ll o, ll m){o%=m;if(o<0)o+=m;return o;}
#define rand my_rand
#define srand my_srand
#define random_shuffle(r...)random_shuffle(r,[](int _){return my_rand()%_;})
sim>inline c nux(c m){if(!m)return numeric_limits<c>::max();c A=m&-m;c B=~((A-1)^m);c C=B&-B;c D=(C>>(1+ctz(A)))-1;return C|(m&~(C-1))|D;}sim>void unq(c&x){x.resize(unique(ALL(x))-x.begin());}
//#STAY AT HOME
//~ https://en.wikipedia.org/wiki/Berlekamp%E2%80%93Massey_algorithm
// const ll mod = (ll)1e18 + 3;
const int mod = 1e9 + 7;
template <class c> void add_self(c & a, c b) { a += b; if(a >= mod) a -= mod; }
template <class c> void sub_self(c & a, c b) { a -= b; if(a < 0) a += mod; }
int mul(int a, int b) { return (ll) a * b % mod; }
ll mul(ll a, ll b) {return (__int128) a * b % mod; }
template <class c> c my_pow(c a, c b) {
	c r = 1;
	while(b) {
		if(b & 1) r = mul(r, a);
		a = mul(a, a);
		b >>= 1;
	}
	return r;
}
template <class c> c my_inv(c a) { return my_pow<c>(a, mod - 2); }
template <class c> c negative (c a) {return (mod - a) % mod;}
template <class F> struct Massey {
	vector<F> start, coef; // 3 optional lines
	// vector<vector<F>> powers;
	F memo_inv;
	
	// Start here and write the next ~25 lines until "STOP"
	
	int L; // L == coef.size() <= start.size()
	Massey(vector<F> in) { // O(N^2)
		L = 0;
		const int N = in.size();
		vector<F> C{1}, B{1};
		for(int n = 0; n < N; ++n) {
			// assert(0 <= in[n] && in[n] < mod); // invalid input
			B.insert(B.begin(), 0);
			F d = 0;
			for(int i = 0; i <= L; ++i)
				add_self(d, mul(C[i], in[n-i]));
			if(d == 0) continue;
			vector<F> T = C;
			C.resize(max(B.size(), C.size()));
			for(int i = 0; i < (int) B.size(); ++i)
				sub_self(C[i], mul(d, B[i]));
			if(2 * L <= n) {
				L = n + 1 - L;
				B = T;
				d = my_inv(d);
				for(F & x : B) x = mul(x, d);
			}
		}
		// cerr << "L = " << L << "\n";
		assert(2 * L <= N - 2); // NO RELATION FOUND :(
		// === STOP ===
		for(int i = 1; i < (int) C.size(); ++i)
			coef.push_back(negative(C[i]));
		assert((int) coef.size() == L);
		for(int i = 0; i < L; ++i)
			start.push_back(in[i]);
		while(!coef.empty() && coef.back() == 0) { coef.pop_back(); --L; }
		if(!coef.empty()) memo_inv = my_inv(coef.back());
		// powers.push_back(coef);
		//~ debug() << imie(coef);
	}
	
	// vector<F> mul_cut(vector<F> a, vector<F> b) {
		// vector<F> r(2 * L - 1);
		// for(int i = 0; i < L; ++i)
			// for(int j = 0; j < L; ++j)
				// add_self(r[i+j], mul(a[i], b[j]));
		// while((int) r.size() > L) {
			// F value = mul(r.back(), memo_inv); // div(r.back(), coef.back());
			// const int X = r.size();
			// add_self(r[X-L-1], value);
			// for(int i = 0; i < L; ++i)
				// sub_self(r[X-L+i], mul(value, coef[i]));
			// assert(r.back() == 0);
			// r.pop_back();
		// }
		// return r;
	// }
	// F get(ll k) { // O(L^2 * log(k))
		// if(k < (int) start.size()) return start[k];
		// if(L == 0) return 0;
		// k -= start.size();
		// vector<F> vec = coef;
		// for(int i = 0; (1LL << i) <= k; ++i) {
			// if(i == (int) powers.size())
				// powers.push_back(mul_cut(powers.back(), powers.back()));
			// if(k & (1LL << i))
				// vec = mul_cut(vec, powers[i]);
		// }
		// F total = 0;
		// for(int i = 0; i < L; ++i)
			// add_self(total, mul(vec[i], start[(int)start.size()-1-i]));
		// return total;
	// }
};
int phi(int x) {
	int ans = 0;
	for (int i = 1; i <= x; ++i) if (gcd(i, x) == 1) ans++;
	return ans;
}
void solve() {
	const int nax = 5040;
	int n;
	cin >> n;
	map <int, int> by_cou;
	REP(_, n) {
		int x;
		cin >> x;
		by_cou[x]++;
	}
	auto compress = [&](const vll &a) {
		assert((int)a.size() == n);
		vll ans;
		int i = 0;
		for (auto [_, c] : by_cou) {
			ll v = 0;
			REP(_, c) v += a[i++];
			ans.push_back(v % mod);
		}
		assert(i == n);
		return ans;
	};
	auto decompress = [&](const vll &a) {
		assert(a.size() == by_cou.size());
		int i = 0;
		vll ans;
		for (auto [_, c] : by_cou) {
			REP(_, c) ans.push_back(a[i]);
			i++;
		}
		assert((int)ans.size() == n);
		return ans;
	};
	vi where(nax, -1);
	vll diag(by_cou.size());
	{
		int i = 0;
		for (auto [v, _] : by_cou) where[v] = i++;
	}
	vector <pair <int, vi> > todo;
	vpii by_cou_vec(ALL(by_cou));
	for (int d = 1; d <= nax; ++d) {
		vi which;
		for (int i = d; i < nax; i += d)
			if (where[i] != -1)
				which.push_back(where[i]);
		if (!which.empty()) {
			int P = phi(d);
			// debug(imie(P), imie(which));
			int actual_count = 0;
			for (int x : which) actual_count += by_cou_vec[x].second;
			if (actual_count > 1) {
				for (int x : which) diag[x] += P * actual_count;
				todo.emplace_back(-P, move(which));
			}
		}
	}
	diag = decompress(diag);
	// debug(imie(diag));
	// debug(imie(todo));
	vector <ll> muls(n - 1);
	for (auto &x : muls) x = 1 + rand() % (mod - 1);
	ll tot_mul = 1;
	for (ll x : muls) (tot_mul *= x) %= mod;
	ll inv = my_inv(tot_mul);
	// cerr << "inv = " << inv << "\n";
	// vector <tuple <int, int, int> > bef, aft;
	int swaps = 0;
	// if (n - 1 >= 2) {
		// REP(_, 2) {
			// REP(id, 100000) {
				// int type = 0;//rand() & 1;
				// int a, b;
				// do {
					// a = id % (n - 1);//rand() % (n - 1);
					// b = rand() % (n - 1);
				// } while (a == b);
				// bef.emplace_back(type ? -1 : rand() % mod, a, b);
				// if (type) swaps++;
			// }
			// swap(bef, aft);
		// }
	// }
	// auto run = [&](const vector<tuple <int, int, int> > which, vector <ll> &b) {
		// for (auto [t, i, j] : which) {
			// if (t == -1) swap(b[i], b[j]);
			// else (b[i] += b[j] * t) %= mod;
		// }
	// };
	auto multiply = [&](vll full_a) {
		// run(bef, full_a);
		full_a.push_back(0);
		assert((int)full_a.size() == n);
		auto a = compress(full_a);
		vll b(by_cou.size());
		for (auto &[d, which] : todo) {
			ll sum = 0;
			for (int c : which) sum += a[c];
			(sum *= d) %= mod;
			for (int c : which) (b[c] += sum) %= mod;
		}
		b = decompress(b);
		assert((int)b.size() == n);
		assert((int)diag.size() == n);
		REP(i, n) (b[i] += diag[i] * full_a[i]) %= mod;
		b.pop_back();
		REP(i, n - 1) b[i] = (b[i] * muls[i]) % mod;
		// run(aft, b);
		return b;
	};
	REP(i, n - 1) {
		vll a(n - 1);
		a[i] = 1;
		// debug(multiply(a));
	}
	vector <ll> powe(n - 1);
	for (auto &x : powe) x = rand() % mod;
	vector <ll> coefs(n - 1);
	for (auto &x : coefs) x = rand() % mod;
	vll mas;
	REP(_, 2 * n + 2) {
		ll prod = 0;
		REP(i, n - 1) (prod += (coefs[i] * powe[i])) %= mod;
		mas.push_back(prod);
		powe = multiply(powe);
	}
	vi mas_int;
	for (auto &x : mas) mas_int.push_back(fix(x, mod));
	// debug(imie(mas));
	Massey<int> m(mas_int);
	// debug(m.coef);
	auto min_poly = move(m.coef);
	ll ans = min_poly.back();
	if ((n + swaps) % 2) ans *= -1;
	(ans *= inv) %= mod;
	// cerr << "rec length = " << min_poly.size() << "\n";
	if ((int)min_poly.size() != n - 1) cout << "0\n";
	else cout << fix(ans, mod) << "\n";
}
int main() {
	ios_base::sync_with_stdio(0);
	cin.tie(0);
	unt seed = TIME; srand(seed);
	// cerr << "seed = " << seed << "\n";
	int t = 1;
	// scanf("%d", &t);
	REP(_, t) solve();
	exit(0);
}
//#STAY AT HOME