BSOJ5660【BZOJ3309】DZY Loves Math

题目

求

\[\large \sum\limits_{i=1}^n\sum\limits_{j=1}^mf(\gcd(i,j)) \]

其中，我们定义 \(f(x)\) 表示：\(x\) 的所有质因子的最大次数。

BSOJ5660【BZOJ3309】DZY Loves Math

分析

还是推一推柿子吧：

\[\large \begin{split} &\ \ \ \ \ \sum\limits_{i=1}^n\sum\limits_{j=1}^m\sum_{d\vert \gcd(i,j)}f(d)[\gcd(i,j)=1]\\ &=\sum_{d=1}^nf(d)\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}[\gcd(i,j)=1]\\ &=\sum_{d=1}^n f(d)\sum_{k=1}^{\lfloor\frac nd\rfloor}\mu(k)\sum_{i=1}^{\lfloor\frac{n}{kd}\rfloor}\sum_{j=1}^{\lfloor\frac m{kd}\rfloor}1\\ &=\sum_{d=1}^n f(d)\sum_{k=1}^{\lfloor\frac nd\rfloor}\mu(k){\lfloor\frac{n}{kd}\rfloor}{\lfloor\frac m{kd}\rfloor} \end{split} \]

然后还是传统艺能，两次除法分块不好做，于是枚举乘积转化成一个整除分块和卷积的形式：

\[\large \sum_{d=1}^n f(d)\sum_{k=1}^{\lfloor\frac nd\rfloor}\mu(k){\lfloor\frac{n}{kd}\rfloor}{\lfloor\frac m{kd}\rfloor}=\sum_{t=1}^n {\lfloor\frac{n}{t}\rfloor}{\lfloor\frac m{t}\rfloor}\sum_{d\vert t}f(d)\mu(\frac td) \]

后面那个东西有个暴力的讨论发现可以直接筛，详见yyb的博客

前面的就是数论分块，结束。

代码

#include<bits/stdc++.h>
using namespace std;
template <typename T>
inline void read(T &x){
	x=0;bool f=false;char ch=getchar();
	while(!isdigit(ch)){f|=ch=='-';ch=getchar();}
	while(isdigit(ch)){x=x*10+(ch^48);ch=getchar();}
	x=f?-x:x;
	return ;
}
template <typename T>
inline void write(T x){
	if(x<0) x=-x,putchar('-');
	if(x>9) write(x/10);
	putchar(x%10^48);
	return ;
}
#define ll long long
#define ull unsigned long long
#define rep(i,x,y) for(int i=(x);i<=(y);i++)
#define dep(i,y,x) for(int i=(y);i>=(x);i--)
const int N=1e7+5,M=2e5+5;
int n,m,prime[N],cnt;
int k,d[N],low[N];
ll Ans,pre[N];
bool vis[N];
void GetPrimes(int v){
	for(int i=2;i<=v;i++){
		if(!vis[i]) prime[++cnt]=i,d[i]=pre[i]=1,low[i]=i;
		for(int j=1;i*prime[j]<v;j++){
			vis[i*prime[j]]=1;
			if(i%prime[j]==0){
				d[i*prime[j]]=d[i]+1;
				low[i*prime[j]]=low[i]*prime[j];
				if (i==low[i]) pre[i*prime[j]]=1;
				else pre[i*prime[j]]=d[i/low[i]]==d[low[i]*prime[j]]?-pre[i/low[i]]:0;
				break;
			}
			d[i*prime[j]]=1;
			low[i*prime[j]]=prime[j];
			pre[i*prime[j]]=d[i]==1?-pre[i]:0;
		}
	}
	for(int i=2;i<=v;i++) pre[i]+=pre[i-1];
	return ;
}
signed main(){
	int t;read(t);GetPrimes(1e7);
	while(t--){
		read(n),read(m);Ans=0;
		if(n>m) swap(n,m);
		for(int l=1,r;l<=n;l=r+1){
			r=min(n/(n/l),m/(m/l));	
			Ans+=1ll*(n/l)*(m/l)*(pre[r]-pre[l-1]);
		}
		write(Ans);putchar('\n');
	}
	return 0;
}