于神之怒加强版

于神之怒加强版

题解

很明显的一道莫比乌斯反演。

原式：\sum_{i=1}^{n}\sum_{j=1}^{m}(i,j)^k

= \sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d=1}^{n}d^k[(i,j)==1]

= \sum_{d=1}^{n}d^k\sum_{j=1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor \frac{m}{d} \rfloor}\sum_{x|i,x|j}\mu\left(x \right )

= \sum_{d=1}^{n}d^k\sum_{x=1}^{\frac{n}{d}}\mu\left(x \right )\lfloor \frac{n}{dx} \rfloor ^2

令T= dx,原式：

=\sum_{T=1}^{n}\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor \sum_{d|T}d^k \mu\left( \frac{T}{d} \right )

令g\left(T \right )= \sum_{d|T}a^k\mu\left(\frac{T}{d} \right )，

反演可得g\left(T \right )=\prod _{i=1}^{t}g\left(P_{i}^{x_{i}} \right )

=\prod _{i=1}^{t}P_{i}^{k(x_{i}-1)}\mu\left(P_{i} \right )+P_{i}^{kx_{i}}\mu\left( 1 \right )

= \prod_{i=1}^{t}P_{i}^{k(x_{i}-1)(P_{i}^{k}-1)}

于是就可以快乐的分块了

源码

#include<cstdio>
#include<cmath>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<queue>
#include<vector>
using namespace std;
#define MAXN 5000010
typedef long long LL;
#define int LL
const int INF=0x7f7f7f7f;
const int mo=1e9+7;
typedef pair<int,int> pii;
#define gc() getchar()
template<typename _T>
void read(_T &x){
	_T f=1;x=0;char s=gc();
	while(s>'9'||s<'0'){if(s=='-')f=-1;s=gc();}
	while(s>='0'&&s<='9'){x=(x<<3)+(x<<1)+(s^48);s=gc();}
	x*=f;
}
int tt,k,n,m,T;
int cntp,prime[MAXN],mobius[MAXN],g[MAXN];
bool oula[MAXN];
int mlt(int a,int b){
	int res=1;
	while(b){
		if(b&1)res=res*a%mo;
		b>>=1;a=a*a%mo;
	}
	return res;
}
void init(){
	mobius[1]=1;
	for(int i=2;i<=5e6;i++){
		if(!oula[i]){
			prime[++cntp]=i;g[cntp]=mlt(i,k);
			mobius[i]=(g[cntp]-1+mo)%mo;
		}
		for(int j=1;j<=cntp&&i*prime[j]<=5e6;j++){
			oula[i*prime[j]]=1;
			if(i%prime[j]==0){
				mobius[i*prime[j]]=mobius[i]*g[j]%mo;
				break;
			}
			mobius[i*prime[j]]=mobius[i]*mobius[prime[j]]%mo;
		}
	}
	for(int i=1;i<=5e6;i++)mobius[i]=(mobius[i-1]+mobius[i])%mo;
}
signed main(){
	read(tt);read(k);init();
	while(tt--){
		read(n);read(m);T=min(n,m);int pos,ans=0;
		for(int i=1;i<=T;i=pos+1){
			pos=min(n/(n/i),m/(m/i));
			ans=(ans+(mobius[pos]-mobius[i-1]+mo)*(n/i)%mo*(m/i)%mo)%mo;
		}
		printf("%lld\n",ans);
	}
	return 0;
}

谢谢！！！

posted @ 2020-03-30 16:06  StaroForgin  阅读(6)  评论(0)    收藏  举报  来源