BZOJ3560 DZY Loves Math V

原题链接：https://www.lydsy.com/JudgeOnline/problem.php?id=3560

DZY Loves Math V

Description

给定n个正整数a1,a2,…,an，求

$$\sum_{i_1|a_1}\sum_{i_2|a_2}...\sum_{i_n|a_n} \varphi(i_1i_2...i_n)$$

的值（答案模10^9+7）。

Input

第一行一个正整数n。

接下来n行，每行一个正整数，分别为a1,a2,…,an。

Output

仅一行答案。

Sample Input

3
6
10
15

Sample Output

1595

HINT

1<=n<=10^5，1<=ai<=10^7。共3组数据。

题解

由于$\varphi$函数为积性函数，不难想到将每个质数的贡献分开计算，最后合并。

对于一个质数$p$，设在$a_i$质因子分解后$p$的指数为$b_i$，那么它的贡献为：

$$\sum_{i_1=0}^{b_1}\sum_{i_2=0}^{b_2}\cdots\sum_{i_n=0}^{b_n}\varphi(p^{\sum_{j=1}^ni_j})$$

因为$\varphi(x)$的定义为小于$x$且与$x$互质的数的个数，那么对于$p^{\sum_{j=1}^ni_j}$来说，与它不互质的数只有$p$的倍数，即每$p$个数就有$p-1$个，所以上式可表示为：

$$\begin{align*} &(\sum_{i_1=0}^{b_1}\sum_{i_2=0}^{b_2}\cdots\sum_{i_n=0}^{b_n}p^{\sum_{j=1}^ni_j}-1)\times \frac{p-1}{p}+1\\ &=(\prod_{i=1}^n\sum_{j=0}^{b_i}p^j-1)\times \frac{p-1}{p}+1 \end{align*}$$

具体实现的话，只需要对每个$a_i$来一发质因数分解，将每个质数在$a_i$中的指数打个包，再以质数为第一关键字，指数为第二关键字排一波序就能方便地处理出上式的前缀和。

代码

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int M=3505,mod=1e9+7;
struct sd{int p,a;}d[M*M];
bool operator <(sd a,sd b){return a.p==b.p?a.a<b.a:a.p<b.p;}
int p[M/3],n,tot;
ll sum[105];
bool isp[M];
void pre(){for(int i=2;i<M;++i){if(!isp[i])p[++p[0]]=i;for(int j=1;j<=p[0]&&i*p[j]<M;++j){isp[i*p[j]]=1;if(i%p[j]==0)break;}}}
void divide(int x)
{for(int i=1;i<=p[0];++i)if(x%p[i]==0){d[++tot]=(sd){p[i],0};for(;x%p[i]==0;x/=p[i],++d[tot].a);}if(x>1)d[++tot]=(sd){x,1};}
ll power(ll x,ll p){ll ans=1;for(;p;p>>=1,x=x*x%mod)if(p&1)ans=ans*x%mod;return ans;}
void in(){int a;pre();scanf("%d",&n);for(int i=1;i<=n;++i)scanf("%d",&a),divide(a);}
void ac()
{
    sort(d+1,d+1+tot);ll ans=1,p,tmp;sum[0]=1;
    for(int i=1,j,k;i<=tot;i=j+1)
    {
        p=d[i].p;tmp=1;
        for(j=i;d[j].p==d[j+1].p&&j<tot;++j);
        for(k=1;k<=d[j].a;++k)sum[k]=sum[k-1]*p%mod;
        for(k=1;k<=d[j].a;++k)(sum[k]+=sum[k-1])%=mod;
        for(k=i;k<=j;++k)tmp=tmp*sum[d[k].a]%mod;
        (ans*=(tmp-1)*(p-1)%mod*power(p,mod-2)%mod+1)%=mod;
    }
    printf("%lld",ans);
}
int main(){in();ac();}