#include <bits/stdc++.h>
using namespace std;
#define ll long long

#define AMOD 1000000007
#define P 599
#define G 7
#define O 17
#define Q (P+2)
#define MOD (P*Q)
#define E ((P*P - 1)/2)
#define EX 19
#define M (1 << EX)
#define MMOD1 13631489LL
#define ROOT1 1306074LL
#define IROOT1 11999696LL
#define SC1 13631463LL
#define MMOD2 14155777LL
#define ROOT2 2742784LL
#define IROOT2 1606624LL
#define SC2 14155750LL
#define IM21 13631463LL
#define IM12 27LL
#define MMOD (MMOD1 * MMOD2)

// arrs
ll *ROOTS1[EX+1];
ll *IROOTS1[EX+1];
ll *ROOTS2[EX+1];
ll *IROOTS2[EX+1];
ll ROOTZ1[EX+1];
ll IROOTZ1[EX+1];
ll ROOTZ2[EX+1];
ll IROOTZ2[EX+1];
int mu[MOD];
bool isprime[MOD];
ll temp[M];
ll ce[M];
ll ci[M];
ll co[M];
ll ge[M];
ll gi[M];
ll go[M];
ll ct[MOD];
ll counts[MOD];
ll gen[E];
ll gin[E];


void fft(ll *a, ll *b, int e, int m, ll **roots, ll mmod) {
    if (e == 0) return;
    int j = 0;
    for (int i = 0; i < m; i += 2) b[j++] = a[i];
    for (int i = 1; i < m; i += 2) b[j++] = a[i];
    ll *w = roots[e--];
    m >>= 1;
    fft(b  , a  , e, m, roots, mmod);
    fft(b+m, a+m, e, m, roots, mmod);
    for (int i = 0, j = m; i < m; i++, j++) {
        ll pd = b[j] * w[i] % mmod;
        a[i] = (b[i] + pd);
        a[j] = (b[i] - pd);
    }
}


void _mulall(ll *ce, ll *ci, ll *co, ll mmod, ll **roots, ll **iroots, ll sc) {
    // compute ce*ce, ce*ci and ci*ci. store in ce, ci and co.

    for (int i = E; i < M; i++) {
        ce[i] = ci[i] = 0;
    }

    fft(ce, temp, EX, M, roots, mmod);
    fft(ci, temp, EX, M, roots, mmod);

    for (int i = 0; i < M; i++) {
        co[i] = ci[i] * ci[i] % mmod;
        ci[i] = ci[i] * ce[i] % mmod;
        ce[i] = ce[i] * ce[i] % mmod;
    }

    fft(ce, temp, EX, M, iroots, mmod);
    fft(ci, temp, EX, M, iroots, mmod);
    fft(co, temp, EX, M, iroots, mmod);

    for (int i = 0; i < M; i++) {
        ce[i] = ce[i] * sc % mmod;
        ci[i] = ci[i] * sc % mmod;
        co[i] = co[i] * sc % mmod;
    }
}


ll chinese(ll a1, ll a2) {
    return (a1 * IM21 % MMOD1 * MMOD2 + a2 * IM12 % MMOD2 * MMOD1) % MMOD;
}


void mulall() {
    // two multiplications done here, one for each of two moduli of the
    // form m*2^k + 1. combine with the chinese remainder theorem.
    for (int i = 0; i < E; i++) {
        ge[i] = ce[i];
        gi[i] = ci[i];
    }
    _mulall(ce, ci, co, MMOD1, ROOTS1, IROOTS1, SC1);
    _mulall(ge, gi, go, MMOD2, ROOTS2, IROOTS2, SC2);
    for (int i = 0; i < M; i++) {
        ce[i] = chinese(ce[i], ge[i]);
        ci[i] = chinese(ci[i], gi[i]);
        co[i] = chinese(co[i], go[i]);
        if (ce[i] < 0) ce[i] += MMOD;
        if (ci[i] < 0) ci[i] += MMOD;
        if (co[i] < 0) co[i] += MMOD;
    }
}

int ob;
void init() {
    // initialize shiz for FFT
    for (int i = 0; i <= EX; i++) {
         ROOTS1[i] = new ll[1 << i];
         ROOTS2[i] = new ll[1 << i];
        IROOTS1[i] = new ll[1 << i];
        IROOTS2[i] = new ll[1 << i];
    }

     ROOTZ1[EX] =  ROOT1;
     ROOTZ2[EX] =  ROOT2;
    IROOTZ1[EX] = IROOT1;
    IROOTZ2[EX] = IROOT2;
    for (int i = EX-1; i >= 0; i--) {
         ROOTZ1[i] =  ROOTZ1[i+1] *  ROOTZ1[i+1] % MMOD1;
        IROOTZ1[i] = IROOTZ1[i+1] * IROOTZ1[i+1] % MMOD1;
         ROOTZ2[i] =  ROOTZ2[i+1] *  ROOTZ2[i+1] % MMOD2;
        IROOTZ2[i] = IROOTZ2[i+1] * IROOTZ2[i+1] % MMOD2;
    }

    for (int i = 0; i <= EX; i++) {
         ROOTS1[i][0] = 1;
         ROOTS2[i][0] = 1;
        IROOTS1[i][0] = 1;
        IROOTS2[i][0] = 1;
        for (int j = 1; j < 1 << i; j++) {
             ROOTS1[i][j] =  ROOTS1[i][j-1] *  ROOTZ1[i] % MMOD1;
             ROOTS2[i][j] =  ROOTS2[i][j-1] *  ROOTZ2[i] % MMOD2;
            IROOTS1[i][j] = IROOTS1[i][j-1] * IROOTZ1[i] % MMOD1;
            IROOTS2[i][j] = IROOTS2[i][j-1] * IROOTZ2[i] % MMOD2;
        }
    }

    int curr = 1, idx = 0;
    do {
        if (curr == O*O) ob = idx;
        gen[idx] = curr;
        gin[idx] = curr * O % MOD;
        idx++;
    } while ((curr = curr * G % MOD) != 1);

    for (int i = 0; i < MOD; i++) {
        mu[i] = 1;
        isprime[i] = true;
    }

    for (int i = 2; i < MOD; i++) {
        if (isprime[i]) {
            for (int j = i; j < MOD; j += i) {
                isprime[j] = false;
                mu[j] = -mu[j];
            }
            if (i <= MOD / i) {
                for (int j = i * i; j < MOD; j += i * i) {
                    mu[j] = 0;
                }
            }
        }
    }
}


// computes the counts array. counts[v] is the number of i,j pairs such that Y(i,j) = v.
void compute_counts() {
    // cases where Y(i,j) == 0 (599) or Y(i,j) == 0 (601) are handled first,
    // as a special case. this one is not hard to optimize into O(MOD) time.

    for (int i = 0; i < MOD; i++) {
        counts[i] = 0;
    }

    ll sp = 0, sq = 0;
    for (int i = 0; i < Q; i++) sp += ct[P * i];
    for (int j = 0; j < P; j++) sq += ct[Q * j];

    for (int i = 0; i < MOD; i++) {
        counts[0] += ct[i] * (1 + (i % P != 0 && i % Q != 0));
    }
    counts[0] -= sp + sq;
    counts[0] *= ct[0];
    counts[0] += sp * sq * 2;

    for (int i = 0; i < Q; i++) {
        temp[i] = 0;
    }
    for (int i = 0; i < MOD; i++) {
        temp[i % Q] += ct[i] * (1 + (i % P != 0 && i % Q != 0));
    }
    for (int i = 1; i < Q; i++) {
        for (int j = P; j < MOD; j += P) {
            counts[i * j % MOD] += temp[i] * ct[j];
        }
    }

    for (int i = 0; i < P; i++) {
        temp[i] = 0;
    }
    for (int i = 0; i < MOD; i++) {
        temp[i % P] += ct[i] * (1 + (i % P != 0 && i % Q != 0));
    }
    for (int i = 1; i < P; i++) {
        for (int j = Q; j < MOD; j += Q) {
            counts[i * j % MOD] += temp[i] * ct[j];
        }
    }

    // starting here is the part that computes the general case
    // O(MOD log MOD) solution using FFT. uses an 'almost-generator'
    // modulo 359999 as base, and does FFT on the exponents.
    // three multiplications are needed; those which can be generated
    // and those not.
    for (int e = 0; e < E; e++) {
        ce[e] = ct[gen[e]];
        ci[e] = ct[gin[e]];
    }

    mulall();

    for (int e = 0; e < E; e++) {
        ce[e] += ce[e+E];
        ci[e] += ci[e+E];
        co[e] += co[e+E];
    }

    for (int e = 0; e < E; e++) {
        counts[gen[e]] += ce[e];
        counts[gin[e]] += ci[e] * 2;
        counts[gen[e]] += co[(e-ob+E)%E];
    }
}


ll solve() {
    compute_counts();

    // given the counts array, compute the triples of coprime number
    // using a 'sieve-like' dynamic programming
    // be careful in this part which mod to use where
    for (int i = 0; i < MOD; i++) {
        counts[i] %= AMOD;
    }

    ll ans = 0;
    ll res = counts[1] * counts[0] % AMOD;
    for (int i = 1; i < MOD; i++) {
        ll dcount = 0;
        for (int j = i; j < MOD; j += i) {
            dcount += counts[j];
        }
        dcount %= AMOD;
        ll pw = dcount;
        res += mu[i] * (pw = pw * dcount % AMOD);
        ans += mu[i] * (pw = pw * dcount % AMOD);
    }

    return (ans + res % AMOD * counts[0] * 3) % AMOD;
}


int main() {
    init();
    int z;
    scanf("%d", &z);
    while (z--) {
        for (int i = 0; i < MOD; i++) {
            ct[i] = 0;
        }
        int n;
        scanf("%d", &n);
        while (n--) {
            int v;
            scanf("%d", &v);
            ct[v % MOD]++;
        }
        ll ans = solve();
        if ((ans %= AMOD) < 0) ans += AMOD;
        printf("%lld\n", ans);
    }
}