spoj GCDEX

传送门：http://www.spoj.com/problems/GCDEX/

思路：令g(n)=ΣΣgcd(i,j)(i<=n,j<=i)，再令f(n)=Σgcd(i,n)

很显然ans=g(n)-(n+1)*n/2

那么g(n)=Σf(i)(i<=n)

f(n)在上一题已证为积性函数，线筛求出，就可以在O(n)时间打出表来。

为了方便先设mindiv[i]为i的最小质因子，T[i]为i的最小质因子次数，d[i]=mindiv[i]^t[i]

具体做法：step1:当n为质数时，显然f(i)=(p-1)*1+1*p=2*p-1

step2:再考虑如何从f(p^a)推出f(p^(a+1))

由上一题的结论得f(p^a)=p^(a-1)*(pa-a+p)

f(p^(a+1))=p*(p^(a-1)*(pa+a-1))

=p*f(p^a)+p^a*(p-1)

然后就是 f[i*pri[j]]=f[i]*pri[j]+d[i]*(pri[j]-1))*f[i/d[i]];

其实如果记了d[i],t[i],f(p^(a+1))有些已经可以直接求（比如这道）....

所以上面这段step2在这题里可不做....

<pre name="code" class="cpp">#include<cstdio>
#include<cstring>
#include<algorithm>
const int N=1000000;
using namespace std;
typedef long long ll;
ll f[N+10],t[N+10],d[N+10],pri[N+10],g[N+10],ans[N+10],tot,x=1;bool isp[N+10];

void prework(){
	memset(isp,1,sizeof(isp));
	isp[1]=0,f[1]=1;
	for (int i=2;i<=N;i++){
		if (isp[i]) pri[++tot]=i,f[i]=i+i-1,d[i]=i,t[i]=1;
		for (int j=1;j<=tot&&pri[j]*i<=N;j++){
			isp[i*pri[j]]=0;
			if (i%pri[j]==0){
				f[i*pri[j]]=f[i]*pri[j]+d[i]*(pri[j]-1))*f[i/d[i]];
				t[i*pri[j]]=t[i]+1,d[i*pri[j]]=d[i]*pri[j];break;
			}
			f[i*pri[j]]=f[i]*f[pri[j]],t[i*pri[j]]=1,d[i*pri[j]]=pri[j];
		}
	}
}

int main(){
	prework();
	g[1]=f[1];for (int i=2;i<=N;i++) g[i]=g[i-1]+f[i];
	ans[1]=1;for (int i=2;i<=N;i++) ans[i]=g[i]-1ll*(i+1)*i/2;
	scanf("%lld",&x);
	while (x) printf("%lld
",ans[x]),scanf("%lld",&x);
	return 0;
}

另一种被0.237s和谐的算法

令g(n)=ΣΣ(gcd(i,j))(i<=n,j<=n)//和上一种不一样

f(n)=(g(n)-(n+1)*n/2)/2

考虑求g(n),有公式Σphi(d)=n (d|n)

g(n)=ΣΣΣphi(d)(1<=i<=n,1<=j<=n,d|gcd(i,j))

即有Σphi(d)*ΣΣ[d|gcd(i,j)](1<=d<=n,1<=i<=n,1<=j<=n)

后面一部分显然是(n/d)^2

于是有Σphi(d)*(n/d)^2

但是回答是O(sqrt(n))的，会被TLE和谐掉