【51Nod 1363】最小公倍数之和(欧拉函数)

题面

传送门

题解

拿到式子的第一步就是推倒

[egin{align} sum_{i=1}^nlcm(n,i) &=sum_{i=1}^nfrac{in}{gcd(i,n)}\ &=nsum_{i=1}^nfrac{i}{gcd(n,i)}\ &=nsum_{d|n}sum_{i=1}^n frac{i}{d}[gcd(n,i)=d]\ &=nsum_{d|n}sum_{i=1}^{frac{n}{d}} i[gcd(i,frac{n}{d})=1]\ end{align} ]

我们定义(varphi(n)=sum_{i=1}^n i[gcd(n,i)=1]),那么上面的式子可以化成(n imes sum_{d|n}varphi(d)),关键在于后半部分

然后来考虑一下这个(varphi)函数有什么性质好了

因为(varphi(n))代表的是所有与(n)互质的数的和……发现和欧拉函数很像诶,那么我们可以猜想一下这个函数说不定也是积性函数?

(sum(n)=sum_{i=1}^ni),考虑对于一个普通的(varphi(n))来说,(1)(n)的和为(frac{n(n+1)}{2}),考虑其中的两个质数(p,q)(p|n,q|n)(varphi(n))的值要减去所有(p)的倍数的值,即减去(p imes sum(frac{n}{p})=frac{n(frac{n}{p}+1)}{2}),同理所有(q)的倍数的值为(frac{n(frac{n}{q}+1)}{2}),但是减去这两个之后又多减了(pq)的倍数的值,要加上(frac{n(frac{n}{pq}+1)}{2})

然后有

[egin{align} frac{n(n+1)}{2}-frac{n(frac{n}{p}+1)}{2}-frac{n(frac{n}{q}+1)}{2}+frac{n(frac{n}{pq}+1)}{2} &=frac{n^2}{2}left(1-frac{1}{p}-1-frac{1}{q}+1-frac{1}{pq} ight)\ &=frac{n^2}{2}left(1-frac{1}{p} ight)left(1-frac{1}{q} ight)\ end{align} ]

然后拓展开来之后,设(p_1,p_2,...)(n)的所有质因子,我们发现$$varphi(n)=frac{n^2}{2}(1-frac{1}{p_1})(1-frac{1}{p_2})...$$

分母里的(2)很麻烦,那么我们就设一个(phi(n)=2varphi(n)),其中(phi(1)=1),那么我们神奇地发现(phi)是一个积性函数!

那么我们只要对(phi)求和,加上(varphi(1))(因为理论上来说(phi(1)=2)但我们算成了(1)),然后再除以(2),就是(varphi)的和了

首先几个边界条件都很好求,有(phi(1)=1,phi(p)=p(p-1),phi(p^c)=p^{2c-1}(p-1))

我们要求(sum_{d|n}phi(d)),考虑这么一个想法,先把(n)给质因数分解,然后枚举一下每一个质因子的次数因为这是个积性函数,所以这个质因子这一部分的贡献是独立的。我们要对它的每一次贡献与其它质因子的贡献相加再乘起来,等价于把它所有能产生的贡献都加起来再与其它的贡献乘起来(就是加法对乘法的分配率),这样就能快速计算了

然后没有然后了

ps:我代码里分解质因数的时候终止条件写成(p[i]leq n)了结果活活(T)死……我给忘了只要(p[i] imes p[i]leq n)就行了……

//minamoto
#include<bits/stdc++.h>
#define R register
#define fp(i,a,b) for(R int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(R int i=a,I=b-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
using namespace std;
char buf[1<<21],*p1=buf,*p2=buf;
inline char getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}
int read(){
    R int res,f=1;R char ch;
    while((ch=getc())>'9'||ch<'0')(ch=='-')&&(f=-1);
    for(res=ch-'0';(ch=getc())>='0'&&ch<='9';res=res*10+ch-'0');
    return res*f;
}
char sr[1<<21],z[20];int C=-1,Z=0;
inline void Ot(){fwrite(sr,1,C+1,stdout),C=-1;}
void print(R int x){
    if(C>1<<20)Ot();if(x<0)sr[++C]='-',x=-x;
    while(z[++Z]=x%10+48,x/=10);
    while(sr[++C]=z[Z],--Z);sr[++C]='
';
}
const int N=1e5+5,P=1e9+7,inv2=500000004;
inline int add(R int x,R int y){return x+y>=P?x+y-P:x+y;}
inline int dec(R int x,R int y){return x-y<0?x-y+P:x-y;}
inline int mul(R int x,R int y){return 1ll*x*y-1ll*x*y/P*P;}
int ksm(R int x,R int y){
	R int res=1;
	for(;y;y>>=1,x=mul(x,x))if(y&1)res=mul(res,x);
	return res;
}
bitset<N>vis;int p[N],v[N],c[N],top,m,sqr,n,res,g,tmp,sum;
void init(int n){
	fp(i,2,n){
		if(!vis[i])p[++m]=i;
		for(R int j=1;j<=m&&1ll*i*p[j]<=n;++j){
			vis[i*p[j]]=1;
			if(i%p[j]==0)break;
		}
	}
}
void solve(int n){
	top=0,res=1,g=n;
	for(R int i=1;i<=m&&1ll*p[i]*p[i]<=g;++i)if(g%p[i]==0){
		v[++top]=p[i],c[top]=0;
		while(g%p[i]==0)g/=p[i],++c[top];
	}
	if(g!=1)v[++top]=g,c[top]=1;
	fp(i,1,top){
		sum=1,g=mul(v[i],v[i]-1);
		for(R int j=1;j<=c[i];++j,g=1ll*g*v[i]%P*v[i]%P)sum=add(sum,g);
		res=mul(res,sum);
	}
	res=mul(res+1,inv2);
	res=mul(res,n);
	print(res);
}
int main(){
//	freopen("testdata.in","r",stdin);
//	freopen("testdata.out","w",stdout);
	int T=read();init(sqr=N-5);
	while(T--)n=read(),solve(n);
	return Ot(),0;
}
原文地址:https://www.cnblogs.com/bztMinamoto/p/10432514.html