51nod1363-最小公倍数之和

Description

给出一个n，求1-n这n个数，同n的最小公倍数的和。
例如：n = 6，1,2,3,4,5,6 同6的最小公倍数分别为6,6,6,12,30,6，加在一起 = 66。
由于结果很大，输出Mod 1000000007的结果。

Input

第1行：一个数T，表示后面用作输入测试的数的数量。（1 <= T <= 50000)
第2 - T + 1行：T个数A[i](A[i] <= 10^9)

Output

共T行，输出对应的最小公倍数之和

Solution

推式子:

[egin{aligned} ans &=sum_{i=1}^{n}{frac{icdot n}{gcd(i,n)}} \ &=ncdot sum_{d|n}{sum_{i=1}^{n/d}i[gcd(i,n/d)==1]} \ &=ncdot sum_{d|n}{sum_{i=1}^{d}i[gcd(i,d)==1]} \ end{aligned} ]

由于(sum_{i=1}^n i cdot [(i,n)=1] = egin{cases} 1 & (n=1) \ frac{phi(n)cdot n}2 & (n>1) end{cases}),

[egin{aligned} ans &=ncdot sum_{d|n}{sum_{i=1}^{d}i[gcd(i,d)==1]} \ &=ncdot (1+frac{sum_{d|n,d>1}{phi(d)cdot d}}{2})\ &=ncdot frac{1+sum_{d|n}{phi(d)cdot d}}2\ end{aligned} ]

考虑 (T(n)=sum_{d|n}phi(d)cdot d), 它是(phi(n) cdot n)和(I)的卷积, 是一个积性函数, 可以通过分解质因数的方式求值.

不妨设(n=prod_{i=1}^B p_i^{a_i}), 再考虑

[egin{aligned} T(p^k) &=1+sum_{i=1}^k p^{k-1}cdot (p-1) cdot p^k \ &=1+sum_{i=1}^k p^{2k-1}cdot (p-1) \ &=1+(p-1)frac{p^{2k+1}-p}{p^2-1} \ &=1+frac{p^{2k+1}-p}{p+1} end{aligned} ]

[T(n)=prod_{i=1}^{B} T(p_i^{a_i}) ]

我们可以 (O(lg^2 n)) 求出(T(n)), 进而求出 (ans).

考虑到质因数分解的复杂度, 用素数定理可求出总复杂度(上界)大约为(O(sqrt n + Tcdotfrac { sqrt n}{lg n}))

Code

#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<queue>
#include<map>
#include<set>
using namespace std;
#define rep(i,l,r) for(register int i=(l);i<=(r);++i)
#define repdo(i,l,r) for(register int i=(l);i>=(r);--i)
#define il inline
typedef long long ll;
typedef double db;

//---------------------------------------
const ll nmod=1000000007,sqnsz=1e5+50;
ll t,n,ans;

int nopr[sqnsz],pr[sqnsz],pp=0;
void init(int bnd){
	nopr[1]=1;
	rep(i,2,bnd){
		if(nopr[i]==0)pr[++pp]=i;
		rep(j,1,pp){
			if(i*pr[j]>bnd)break;
			nopr[i*pr[j]]=1;
			if(i%pr[j]==0)break;
		}
	}
}
int fac[sqnsz][2],pf=0;
void getfac(ll n){
	int tmp=sqrt(n),tmp1;
	rep(i,1,pp){
		tmp1=pr[i];
		if((ll)tmp1*tmp1>n)break;
		if(n%tmp1==0){
			fac[++pf][0]=tmp1,fac[pf][1]=0;
			while(n%tmp1==0)++fac[pf][1],n/=tmp1; 
		}
	}
	if(n!=1)fac[++pf][0]=n,fac[pf][1]=1;
}

ll qp(ll a,ll b){
	ll res=1;
	for(;b;a=a*a%nmod,b>>=1){
		if(b&1)res=res*a%nmod;
	}
	return res;
}

int main(){
	ios::sync_with_stdio(0),cin.tie(0);
	init(1e5+5);
	cin>>t;
	rep(cs,1,t){
		cin>>n;
		pf=0,getfac(n);
		ans=1;
		ll a,b,tmp;
		rep(i,1,pf){
			a=fac[i][0],b=fac[i][1];
			tmp=(1+(qp(a,b*2+1)-a)*qp(a+1,nmod-2))%nmod;
			ans=ans*tmp%nmod;
		}
		ans=(ans+1)*n%nmod*qp(2,nmod-2)%nmod;
		if(ans<0)ans+=nmod;
		cout<<ans<<'
';
	}
	return 0;
}