#include <sstream>
#include <string>
#include <vector>
#include <map>
#include <algorithm>
#include <iostream>
#include <utility>
#include <set>
#include <cctype>
#include <queue>
#include <stack>
#include <cstdio>
#include <cstdlib>
#include <cmath>

#define INF 1000000000

using namespace std;
typedef long long ll;
const int MAXK = 100100;
const int MAXT = 1000010;
const ll MOD = 1000000007;
int K;
int C[MAXK];
ll fact[MAXT], inv[MAXT], normalinv[MAXT];

ll modpow(ll a, ll n) {
    if (n == 0) return 1;
    if (n == 1) return a;
    if (n % 2 == 0) {
        ll tmp = modpow(a, n/2);
        return (tmp * tmp) % MOD;
    } else {
        ll tmp = modpow(a, n-1);
        return (tmp * a) % MOD;
    }
}

// extgcd
ll extgcd(ll a, ll b, ll& x, ll& y) {
    ll d = a;
    if (b != 0) {
        d = extgcd(b, a % b, y, x);
        y -= (a / b) * x;
    } else {
        x = 1; y = 0;
    }
    return d;
}

// mod_inverse
ll mod_inverse(ll a, ll m = MOD) {
    ll x, y;
    extgcd(a, m, x, y);
    return (m+x%m) % m;
}

void init() {
    fact[0] = inv[0] = 1;
    for (ll i = 1; i < MAXT; i++) {
        fact[i] = (fact[i-1] * i) % MOD;
        inv[i] = mod_inverse(fact[i]);
    }
    for (ll i = 1; i < MAXT; i++) {
        normalinv[i] = mod_inverse(i);
    }
}

int main(void) {
    init();
    //for (int i = 0; i < 10; i++) {
    //    cout << fact[i] << "  " << inv[i] << endl;
    //}
    cin >> K;
    int g = 0, T = 0;
    for (int i = 0; i < K; i++) {
        scanf("%d", C+i);
        g = __gcd(g, C[i]);
        T += C[i];
    }
    vector<int> V;
    for (int i = 1; i * i <= g; i++) {
        if (g % i == 0) {
            V.push_back(T/i);
            if (i*i != g) V.push_back(T/(g/i));
        }
    }
    sort(V.begin(), V.end());
    ll ans = 0;
    vector<ll> P(V.size());
    for (int i = 0; i < V.size(); i++) {
        ll G = T / V[i];
        P[i] = fact[V[i]];
        for (int j = 0; j < K; j++) {
            P[i] = (P[i] * inv[C[j]/G]) % MOD;
        }
        for (int j = 0; j < V.size(); j++) {
            if (j >= i) break;
            if (V[i] % V[j] == 0) P[i] = (P[i]-P[j]) % MOD;
            if (P[i] < 0) P[i] += MOD;
        }
        ans = (ans + P[i] * normalinv[V[i]]) % MOD;
    }
    cout << ans << endl;
    return 0;
}