[bzoj3560] DZY Loves Math V

Description###

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

[sumlimits_{i_1|a_1} sumlimits_{i_2|a_2} … sumlimits_{i_n|a_n} varphi(i_1 i_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(i)) 为积性函数
设 (i_1 i_2 ··· i_n = p_1^{r_1} p_2^{r_2} … p_m^{r_m}) (p为质数)，那么 (varphi(i_1 i_2 ··· i_n)=varphi(p_1^{r_1}) varphi(p_2^{r_2}) … varphi(p_m^{r_m}))
那么可以对于每个p单独考虑它对答案的贡献，最后都乘起来就行了

乱入
在(i_1,i_2,···,i_n)互不干扰时

[sumlimits_{i_1} sumlimits_{i_2} … sumlimits_{i_n} i_1 i_2 … i_n = (sum_{i_1})(sum_{i_2}) … (sum_{i_n}) ]

上面所说的“乘起来”就是这个道理。
于是对于某个p，设它在(a_i)中的指数为(b_i)。那它对于答案的贡献为

[egin{equation*} egin{aligned} &sumlimits_{i_1=0}^{b_1} sumlimits_{i_2=0}^{b_2}… sumlimits_{i_n=0}^{b_n} varphi(p^{i_1+i_2+…+i_n}) \ =&frac{p-1}{p}[(sumlimits_{i_1=0}^{b_1} sumlimits_{i_2=0}^{b_2} … sumlimits_{i_n=0}^{b_n} p^{i_1+i_2+…+i_n})-1]+1 \ =&frac{p-1}{p}[(sumlimits_{i_1=0}^{b_1} sumlimits_{i_2=0}^{b_2} … sumlimits_{i_n=0}^{b_n} p^{i_1}p^{i_2}…p^{i_n})-1]+1 \ =&frac{p-1}{p}[(p^0+p^1+…+p^{b_1})(p^0+p^1+…+p^{b_2})…(p^0+p^1+…+p^{b_n})-1]+1 end{aligned} end{equation*} ]

(第一个等号后“+1 -1”是考虑1的情况；第三个等号用了上面乱入的道理)

然后就比较好办了。
将 (a_i) 分解质因数，对于它的每个质因子p ，累计p对答案的贡献
最后统一计算

---

代码##

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>

#define P 1000000007

using namespace std;

typedef long long ll;
const int N = 100005;
const int M = 10000005;

int pre[M],prime[M],pnum;
void getp(){
    pre[1]=1;
    for(int i=2;i<M;i++){
        if(!pre[i]) pre[i]=prime[pnum++]=i;
        for(int j=0;j<pnum && (ll)i*prime[j]<M;j++){
            pre[i*prime[j]]=prime[j];
			if(i%prime[j]==0) break;
        }
    } 
}

int n,top;
int a[N],st[N],cnt[N];
ll p[M];

ll Pow_mod(ll x,ll y){
	ll ret=1;
	while(y){
		if(y&1) ret=(ret*x)%P;
		x=(x*x)%P;
		y>>=1;
	}
	return ret;
}
ll inv(ll x) { return Pow_mod(x,P-2); }

int main()
{
 	ll sum,last;
    getp();
    for(int i=0;i<M;i++) p[i]=1;
    scanf("%d",&n);
    for(int i=0;i<n;i++){
        scanf("%d",&a[i]);
		top=0;
		for(int j=a[i];j>1;j=j/pre[j]){ //分解质因数（新技能qwq）
		    if(pre[j]!=st[top]) st[++top]=pre[j];
			cnt[top]++;
		}
		for(int j=1;j<=top;j++){
			sum=1; last=1;
		    while(cnt[j]--) { sum=(sum+last*st[j])%P; last=(last*st[j])%P; }		
			p[st[j]]=(p[st[j]]*sum)%P;
			cnt[j]=0;
		}
    }
    
    ll ans=1;
    for(int i=2;i<M;i++){
		if(p[i]==1) continue;
		p[i]=(((p[i]-1)*(i-1)%P)*inv((ll)i)%P+1)%P;
		ans=(ans*p[i])%P;
	}
	printf("%lld
",ans);
    
    return 0;
}