[BZOJ2693]jzptab

Description
\(\sum\limits_{i=1}^n\sum\limits_{j=1}^m lcm(i,j)\),答案模1e9+9输出,多组询问

Input
一个正整数T表示数据组数
接下来T行 每行两个正整数 表示N、M

Output
T行 每行一个整数 表示第i组数据的结果

Sample Input
1
4 5

Sample Output
122

HINT
T <= 10000
N, M<=10000000


我们令n<m,然后将柿子化简

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m lcm(i,j) \]

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m \dfrac{i\times j}{\gcd(i,j)} \]

\[\sum\limits_{d=1}^n\sum\limits_{i=1}^n\sum\limits_{j=1}^m\dfrac{i\times j}{d}[\gcd(i,j)=d] \]

\[\sum\limits_{d=1}^n\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}\dfrac{d^2\times i\times j}{d}[\gcd(i,j)=1] \]

\[\sum\limits_{d=1}^n d\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}i\times j\sum\limits_{x|i,x|j}\mu(x) \]

\[\sum\limits_{d=1}^n d\sum\limits_{x=1}^{\lfloor\frac{n}{d}\rfloor}\mu(x)\times x^2\sum\limits_{i=1}^{\lfloor\frac{n}{dx}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{dx}\rfloor}i\times j \]

我们发现最后那个是等差数列,继续化简

\[\sum\limits_{d=1}^n d\sum\limits_{x=1}^{\lfloor\frac{n}{d}\rfloor}\mu(x)\times x^2\dfrac{\lfloor\frac{n}{dx}\rfloor(\lfloor\frac{n}{dx}\rfloor+1)\lfloor\frac{m}{dx}\rfloor(\lfloor\frac{m}{dx}\rfloor+1)}{4} \]

然后我们令\(T=dx\),那么得到

\[\sum\limits_{T=1}^n\dfrac{\lfloor\frac{n}{T}\rfloor(\lfloor\frac{n}{T}\rfloor+1)\lfloor\frac{m}{T}\rfloor(\lfloor\frac{m}{T}\rfloor+1)}{4}T\sum\limits_{x|T}\mu(x)x \]

我们设\(f(T)=\sum\limits_{x|T}\mu(x)x\),预处理出f,就可以分块了

\(g(x)=\mu(x)x\),当a,b互质,\(g(a)\times g(b)=ab\mu(a)\mu(b)=ab\mu(ab)=g(ab)\),所以g是积性函数,根据莫比乌斯反演的性质,f也是积性函数

\(T=\prod\limits_{i=1}^k P_i^{x_i}\)\(f(P_i^{x_i})=(1-P_i)\),那么

\[f(T)=\prod\limits_{i=1}^k(1-P_i) \]

这样子我们可以在\(O(n)\)时间内线筛出来,然后维护一下\(T\times f(T)\)的前缀和,然后就可以在\(O(\sqrt N)\)的时间内完成每次询问

/*program from Wolfycz*/
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define inf 0x7f7f7f7f
using namespace std;
typedef long long ll;
typedef unsigned int ui;
typedef unsigned long long ull;
inline int read(){
	int x=0,f=1;char ch=getchar();
	for (;ch<'0'||ch>'9';ch=getchar())	if (ch=='-')    f=-1;
	for (;ch>='0'&&ch<='9';ch=getchar())	x=(x<<1)+(x<<3)+ch-'0';
	return x*f;
}
inline void print(int x){
	if (x>=10)	print(x/10);
	putchar(x%10+'0');
}
const int N=1e7,p=1e8+9;
int prime[N+10],f[N+10],tot;
bool inprime[N+10];
void prepare(){
	f[1]=1;
	for (int i=2;i<=N;i++){
		if (!inprime[i])	prime[++tot]=i,f[i]=1-i+p;
		for (int j=1;j<=tot&&i*prime[j]<=N;j++){
			inprime[i*prime[j]]=1;
			if (i%prime[j]==0){
				f[i*prime[j]]=f[i];
				break;
			}
			f[i*prime[j]]=1ll*f[i]*f[prime[j]]%p;
		}
	}
	for (int i=1;i<=N;i++)	f[i]=(f[i-1]+1ll*i*f[i]%p)%p;
}
int get(int x){return (1ll*x*(x+1)>>1)%p;}
int main(){
	prepare();
	for (int Data=read();Data;Data--){
		int n=read(),m=read(),pos,Ans=0;
		if (n>m)	swap(n,m);
		for (int T=1;T<=n;T=pos+1){
			pos=min(n/(n/T),m/(m/T));
			Ans=(Ans+1ll*get(n/T)*get(m/T)%p*(f[pos]-f[T-1]+p)%p)%p;
		}
		printf("%d\n",Ans);
	}
	return 0;
}
作者:Wolfycz
本文版权归作者和博客园共有,欢迎转载,但必须在文章开头注明原文出处,否则保留追究法律责任的权利
posted @ 2018-08-16 20:59  Wolfycz  阅读(143)  评论(0编辑  收藏  举报