【BZOJ2693】jzptab

题目大意

给定(n,m,(1leq n,m leq 1e7))，求(sumlimits_{i=1}^nsumlimits_{j=1}^mlcm(i,j))，答案对(1e8+9)取模，(T(1leq T leq 1e4))组询问。

题目分析

前半部分同【2011集训贾志鹏】Crash 的数字表格。

但是，由于多组询问，原来的(O(n))做法已经无法满足要求，所以，我们需要更优秀的(O(sqrt n))做法。

以下式子从【2011集训贾志鹏】Crash 的数字表格的最后一步开始化简。

[egin{split} ans &=sumlimits_{d=1}^n dcdot f(1)\ &=sumlimits_{d=1}^n dcdot sumlimits_{i=1}^{lfloorfrac n d floor}mu(i)cdot g(i)\ &=sumlimits_{d=1}^n dcdot sumlimits_{i=1}^{lfloorfrac n d floor}mu(i)cdot icdot icdot sum(lfloorfrac n {di} floor)cdot sum(lfloorfrac m {di} floor)\ &=sumlimits_{i=1}^nmu(i)cdot i^2sumlimits_{d=1}^{lfloorfrac ni floor}sum(lfloorfrac n{di} floor)cdot sum(lfloorfrac m{di} floor)cdot d\ &=sumlimits_{i=1}^nmu(i)cdot i^2sumlimits_{i|T}^nsum(lfloorfrac nT floor)cdot sum(lfloorfrac mT floor)cdot frac Ti\ &=sumlimits_{T=1}^nsum(lfloorfrac nT floor)cdot sum(lfloorfrac mT floor)sumlimits_{i|T}mu(i)cdot i^2cdot frac Ti\ &=sumlimits_{T=1}^nsum(lfloorfrac nT floor)cdot sum(lfloorfrac mT floor)sumlimits_{i|T}mu(i)cdot icdot T\ &=sumlimits_{T=1}^nsum(lfloorfrac nT floor)cdot sum(lfloorfrac mT floor)cdot Tsumlimits_{i|T}mu(i)cdot i\ end{split} ]

其中(sum(i)=frac{(1+i)cdot(i)}{2}),可以直接计算出；(sumlimits_{i|T}mu(i)cdot i)不是积性函数，但仍可以预处理出。

套上整除分块，这个算法的时间复杂度便达到了(O(sqrt n))，可以满足题目要求。

代码实现

#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdio>
#include<iomanip>
#include<cstdlib>
#define MAXN 0x7fffffff
typedef long long LL;
const int N=1e7+5,mod=1e8+9;
using namespace std;
inline int Getint(){register int x=0,f=1;register char ch=getchar();while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}return x*f;}
int g[N],prime[N];
bool vis[N];
int sum(int x){return 1ll*(1+x)*x/2%mod;}
int main(){
	g[1]=1;
	for(int i=2;i<=1e7;i++){
		if(!vis[i])prime[++prime[0]]=i,g[i]=1-i;
		for(int j=1;j<=prime[0]&&1ll*i*prime[j]<=1e7;j++){
			vis[i*prime[j]]=1;
			if(i%prime[j]==0){
				g[i*prime[j]]=g[i];
				break;
			}
			g[i*prime[j]]=(g[i]-1ll*prime[j]*g[i])%mod;
		}
	} 
	for(int i=1;i<=1e7;i++)g[i]=(1ll*g[i]*i+g[i-1])%mod;
	
	int T=Getint();
	while(T--){
		int n=Getint(),m=Getint();
		if(n>m)swap(n,m);
		int ans=0;
		for(int l=1,r;l<=n;l=r+1){
			r=min(n/(n/l),m/(m/l));
			ans=(ans+1ll*sum(n/l)*sum(m/l)%mod*(g[r]-g[l-1])%mod)%mod; 
		}
		cout<<(ans+mod)%mod<<'
';
	}
	return 0;
}