1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
#include <bits/stdc++.h>
using namespace std;
#define rep(i,a,n) for (int i=a;i<n;i++)
#define per(i,a,n) for (int i=n-1;i>=a;i--)
#define pb push_back
#define eb emplace_back
#define mp make_pair
#define all(x) (x).begin(),(x).end()
#define fi first
#define se second
#define SZ(x) ((int)(x).size())
typedef vector<int> VI;
typedef basic_string<int> BI;
typedef long long ll;
typedef pair<int,int> PII;
typedef double db;
mt19937 mrand(1);
const ll mod=1000000007;
int rnd(int x) { return mrand() % x;}
ll powmod(ll a,ll b) {ll res=1;a%=mod; assert(b>=0); for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;}
//ll gcd(ll a,ll b) { return b?gcd(b,a%b):a;}
// head

const int M=5000;
const int N=5010;
int n,a[N],phi[N],cc[N],c[N],pos[N];
VI d[N];
ll f[N][N],ff[N];

ll det(int n) {
	ll ans=1;
	per(i,1,n+1) {
		if (f[i][i]==0) {
			per(j,1,i) if (f[j][i]!=0) {
				rep(k,1,i+1) swap(f[i][k],f[j][k]);
				ans*=-1;
				break;
			}
		}
		if (f[i][i]==0) return 0;
		ll v=powmod(mod-f[i][i],mod-2);
		int t=0;
		rep(k,1,i+1) if (f[i][k]) {
			pos[t++]=k;
		}
		rep(j,1,i) if (f[j][i]) {
			ll tmp=f[j][i]*v%mod;
			rep(kk,0,t) {
				int k=pos[kk];
				f[j][k]=(f[j][k]+tmp*f[i][k])%mod;
			}
		}
	}
	rep(i,1,n+1) ans=ans*f[i][i]%mod;
	if (ans<0) ans+=mod;
	return ans;
}
int main() {
	scanf("%d",&n);
	//n=5000;
	rep(i,0,n) {
		scanf("%d",&a[i]);
		//a[i]=i+1;
	}
	for (int i=1;i<=M;i++) {
		phi[i]+=i;
		for (int j=2*i;j<=M;j+=i) phi[j]-=phi[i];
	}
	for (int i=1;i<=M;i++) for (int j=i;j<=M;j+=i)
		d[j].pb(i);
	rep(i,0,n) for (auto v:d[a[i]]) cc[v]+=1;
	rep(i,1,n) {
		for (auto v:d[a[i]]) c[i]+=cc[v]*phi[v];
		//printf("%d %d\n",i,c[i]);
		ll inv=powmod(c[i],mod-2);
		for (auto v:d[a[i]]) ff[v]=(ff[v]+inv)%mod;
	}
	ll ans=1;
	rep(i,1,n) ans=ans*c[i]%mod;
	rep(i,1,M+1) {
		f[i][i]+=1;
		for (int dd:d[i]) for (int x=1;x*i<=M;x+=1)
			if (gcd(x,i/dd)==1) {
				int j=x*dd;
				f[i][j]=(f[i][j]-phi[i]*ff[lcm(i,j)])%mod;
			}
	}
	//printf("%d\n",(int)clock());
	ll v=det(M);
	//printf("%lld\n",v);
	ans=ans*v%mod;
	if (ans<0) ans+=mod;
	printf("%lld\n",ans);
}