于神之怒加强版

2024-03-16 23:38

文章标签 加强版神之怒

本文主要是介绍于神之怒加强版，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

于神之怒加强版

题解

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

原式： $\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;
}