容易注意到每个质因数是独立的，求出来加起来即可。那就分解质因数对每个质因数搞一个序列，每个值为对应数的该质因数次数，求答案。这样非零数的总数是线对，但是加上零就爆炸了，所以我们的复杂度要严格只与非零数的个数相关。另：分解质因数并不需要线根的做法，普通分解质因数的结果容量是对数，复杂度瓶颈在于每次找最小因数，这是根号的。我们可以预先线对地预处理每个数的最小非一因数就能做到总线对了。
现在考虑对某个质数如何求答案。其实就是选一个最终到达的点 (x)，使 (sum|v_i-x|) 最小。那么跟据某些结论，(x) 肯定取 (v) 的中位数。但是我们不关心中位数，只关心上述和式的值，那显然就是左一半是 (-v_i)，右一半 (v_i)，中间（如果有的话）是 (0)。我们要对 (2^n) 个子集（包括那些占了大量空间的零）求其和。于是我们考虑计算每个数的贡献，也就是 (v_i) 乘上在右一半的次数减去在左一半的次数。这时候 (v_0=0) 显然没有贡献，这就（一定程度上）保证了复杂度。
那考虑排序（可以桶排（因值很小）来维持总 1log），如果 (v_x) 左边选的数严格小于右边选的数，那么就属于左半边；右半边类似。现在来求出现在左半边的次数，显然就是 (sumlimits_{i=0}^{x-1}sumlimits_{j=i+1}^{n-x}dbinom{x-1}idbinom{n-x}{j})。如果对每个位置都平方求其值，那么复杂度是三方对数的。我们考虑 (mathrm O(1)) 求。
注意到 ((x-1)+(n-x)=n-1) 是常数，不难想到范德蒙德卷积。那么我们要将 (i+j=k) 分出每个 (k) 相等的类，并且最好能填满，我们来试试行不行。将可行的 (i,j) 画到二维表上发现是个主对角线的上三角，那就可以分出若干 (i-j) 相等的满类，这不是我们想要的。但注意到 (dbinom ab=dbinom a{a-b})，可以转差为和。
下一步枚举 (i+j=k) 的 (k)。显然合法的 (kin[0,i+n-x-i-1=n-x-1])。然后枚举 (i)，要满足 (iin[0,min(k,x-1)]) 且 (jin[0,n-x-i-1]) 即 (k-iin[0,n-x-i-1])，容易发现 (k-ileq n-x-i-1) 就跟没说一样，剩下 (k-igeq 0) 即 (ileq k)。总结起来就是 (iin[0,min(k,x-1)])。。。。
容易发现 (i) 已经将 ([0,min(k,x-1,n-x)]) 填满了（越界的都是 (0)），那么后面的和式可以直接范德蒙德卷积。所以 (ans=sumlimits_{i=0}^{n-x-1}dbinom{n-1}i)，这就预处理一个第 (n-1) 行的组合数前缀和就可以了，(mathrm O(1)) 目标达成。关于右半边，用总数减去左半边和恰好中间就好了。中间就是 (sumlimits_idbinom{x-1}idbinom{n-x}i=sumlimits_idbinom{x-1}idbinom{n-x}{n-x-i}=dbinom{n-1}{n-x})，易得右半边是 (sumlimits_{i=n-x+1}^{n-1}dbinom{n-1}i)。于是这题就做完啦！总复杂度线对。
code
#include<bits/stdc++.h>
using namespace std;
#define pb push_back
const int N=1000010,inf=0x3f3f3f3f;
const int mod=1e9+7;
int qpow(int x,int y){
int res=1;
while(y){
if(y&1)res=1ll*res*x%mod;
x=1ll*x*x%mod;
y>>=1;
}
return res;
}
int inv(int x){return qpow(x,mod-2);}
int n;
int least[N];
vector<int> buc[N];
int fact[N],factinv[N];
int comb(int x,int y){return 1ll*fact[x]*factinv[y]%mod*factinv[x-y]%mod;}
int C[N];
int cnt[30];
int main(){
cin>>n;
fact[0]=factinv[0]=1;for(int i=1;i<=n;i++)factinv[i]=inv(fact[i]=1ll*fact[i-1]*i%mod);
C[0]=comb(n-1,0);for(int i=1;i<n;i++)C[i]=(C[i-1]+comb(n-1,i))%mod;
memset(least,0x3f,sizeof(least));
for(int i=2;i<=1e6;i++)for(int j=2*i;j<=1e6;j+=i)least[j]=min(least[j],i);
for(int i=1;i<=n;i++){
int x;
scanf("%d",&x);
while(true){
int cnt=0,j=least[x];
if(j==inf)break;
while(x%j==0)x/=j,cnt++;
buc[j].pb(cnt);
}
if(x>1)buc[x].pb(1);
}
int ans=0;
for(int i=1;i<=1e6;i++)if(buc[i].size()){
vector<int> &v=buc[i];
for(int j=0;j<v.size();j++)cnt[v[j]]++;
int now=0;
for(int j=1;j<30;j++)while(cnt[j])v[now++]=j,cnt[j]--;
for(int j=0;j<v.size();j++){
int x=j+1+(n-v.size());
(ans+=(-1ll*C[n-x-1]*v[j]%mod+1ll*(C[n-1]-C[n-x])*v[j]%mod)%mod)%=mod;
}
}
cout<<(ans+mod)%mod;
return 0;
}