from math import ceil

def gcd(a, b):
    if a > b:
        return gcd(b, a)
    while a != 0:
        a, b = b % a, a
    return b

def product(lst):
    prod = 1
    for n in lst:
        prod = (prod * n) % (10 ** 9 + 7)
    return prod
    
# Have a global variable denote the bounds
# Have the function denote how much to divide the bounds by
bounds = None
cache = None
    
def coprime_count(divisor):
    if divisor > bounds[-1] // 2:
        return 1
    if divisor in cache:
        return cache[divisor]
    # bounds are assumed to be sorted
    scaled_bounds = [n // divisor for n in bounds]
    total = product(scaled_bounds)
    cap = scaled_bounds[0]
    # total = sum d from 1 to lowest bound coprime_count(d) (inclusive)
    rest = sum(coprime_count(divisor * d) for d in range(2, cap + 1))
    cache[divisor] = (total - rest) % (10 ** 9 + 7)
    return cache[divisor]
    
t = int(input())
for _ in range(t):
    k = int(input())
    bounds = list(map(int, input().split()))
    bounds.sort()
    cache = dict()
    total = 0
    for i in range(1, bounds[0]+1):
        total = (total + i * coprime_count(i)) % (10 ** 9 + 7)
    print(total)

//https://www.codechef.com/viewsolution/44723407

#include <bits/stdc++.h>
#define FAST ios_base::sync_with_stdio(false);cin.tie(NULL);cout.tie(NULL);
typedef long long ll;
typedef long double ld;
#define pb push_back
#define mp make_pair
#define ff first
#define ss second
#define mod 1000000007
#define pii pair<ll,ll>
#define inf 1000000000000000000
#define bpc(x) __builtin_popcountll(x)
#define autoit(x,it) for(auto it = x.begin(); it != x.end(); it++)
#define autoitr(x,it) for(auto it = x.rbegin(); it != x.rend(); it++)
#define rep(n) for(ll i = 0; i < n; i++)
#define repi(i,n) for(ll i = 0; i < n; i++)
#define hmap gp_hash_table<ll, ll>

#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
using namespace __gnu_pbds;

#define ordered_set tree<ll, null_type,less<ll>, rb_tree_tag,tree_order_statistics_node_update>

using namespace std;
mt19937_64 mt(chrono::steady_clock::now().time_since_epoch().count());

vector<ll> mu, phi, primes, f;

void calc(ll n)
{
    mu.resize(n + 1);
    phi.resize(n + 1);
    f.resize(n + 1);
    primes.clear();
    vector<ll> nums(n + 1);
    nums[1] = 1;
    f[1] = 1;

    bool isp[n + 1];
    rep(n + 1)
    isp[i] = true;
    isp[0] = isp[1] = false;
    mu[1] = 1, phi[1] = 1;
    for (ll i = 2; i <= n; i++)
    {
        if (isp[i])
        {
            phi[i] = i - 1;
            primes.pb(i);
            mu[i] = -1;
            nums[i] = i;
            f[i] = (i - 1)%mod;               ////!!!!!HANDLE THIS CASE FOR A FUNCTION
            //When number is a prime
        }

        ll tmp = n / i;

        for (ll j = 0; j < primes.size() && primes[j] <= tmp; j++)
        {
            ll curr = primes[j] * i;
            isp[curr] = false;
            if (i % primes[j] == 0)
            {
                nums[curr] = nums[i] * primes[j];
                phi[curr] = primes[j] * phi[i];
                mu[curr] = 0;
                if (curr == nums[curr])
                    f[curr] = (curr - i)%mod;                           //!!!!!HANDLE THIS CASE FOR A FUNCTION
                //when number is some power(>=2) of primes[j]
                else f[curr] = (f[curr / nums[curr]] * f[nums[curr]])%mod;            //when there are at least 2 primes and power of primes[j] is atleast 2 in curr
                break;
            }

            nums[curr] = primes[j];
            mu[curr] = (mu[i] * mu[primes[j]]);
            phi[curr] = (phi[i] * phi[primes[j]]);
            f[curr] = (f[primes[j]] * f[i])%mod;                                //When number has only single power of this prime
        }
    }

}

ll powa(ll a, ll b, ll c)
{
    a%=c;
    if(a<0)
    a+=c;
    ll res = 1;
    while(b>0)
    {
        if(b&1)
            res*=a, res%=c;
        a*=a, a%=c;
        b>>=1;
    }
    return res;
}
 
ll norm(ll a)
{
    a%=mod;
    if(a<0)
        a+=mod;
    return a;    
}
 
ll add(ll a, ll b)
{
    a = norm(a);
    b = norm(b);
    return norm(a+b);
}
 
ll mul(ll a, ll b)
{
    a = norm(a);
    b = norm(b);
    return norm(a*b);
}
 
 
ll sub(ll a, ll b)
{
    a = norm(a);
    b = norm(b);
    return norm(a-b);
}
 
ll divi(ll a, ll b)
{
    b = powa(b,mod-2,mod);
    return mul(a,b);
}

#define N 100005

int main()
{
    FAST/**/
    
    ll inv[N];
    inv[0] = 0;
    inv[1] = 1;
    for(ll i=2;i<N;i++)
        inv[i] = divi(1, i);

    ll t;
    cin >> t;

    calc(N-1);

    while (t--)
    {
        ll k;
        cin>>k;
        
        ll arr[k];
        ll mx = N;
        rep(k)
                cin>>arr[i], mx = min(mx, arr[i]);
        ll ppr[N];
        rep(N)
                ppr[i] = 1;
        for(ll i=0;i<k;i++)
        {
            ll val = arr[i];
            ll id = 1;
            for(id = 1; id*id <= val; id++)
            {
                ppr[id] = mul(ppr[id], val/id);
                ppr[id+1] = mul(ppr[id+1], inv[val/id]);
            }
            
            //assert(id>0);
            
            for(ll cid = val/id; cid > 0; cid--)
            {
                ll l = (val+1+cid)/(cid+1), r = (val/cid);
                l = max({l, id});
                r = min({r, val});
                if(l>r)
                    continue;
                ppr[l] = mul(ppr[l], cid);
                ppr[r+1] = mul(ppr[r+1], inv[cid]);
            }
        }
        
        for(ll i=2;i<=mx;i++)
            ppr[i] = mul(ppr[i-1], ppr[i]);
        
        ll ans = 0;
        for(ll i=1;i<=mx;i++)
        {
            ll tmp = ppr[i];
            tmp = mul(tmp, f[i]);
            ans = add(ans, tmp);
        }
        cout<<ans<<"\n";
                

    }

    return 0;
}