bzoj4407:於神之怒加強版
傳送門
這個題真的也是有點難度啊（應該是因為我太菜了）
\[
ans=\sum_{i=1}^{n}\sum_{j=1}^{m}gcd(i,j)^k\\]
可以設
\[
f(d)=\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)==d]d^k\ans=\sum_{d=1}^{min(n,m)}f(d)=\sum_{d=1}^{min(n,m)}\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)==d]d^k\ans=\sum_{d=1}^{min(n,m)}d^k\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)==d]\\]
設
\[ F(d)=\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)==d]\g(x)=\sum_{x|d}F(d)=\sum_{i=1}^{n}\sum_{j=1}^{m}[x|gcd(i,j)]=\lfloor \frac{n}{x}\rfloor\lfloor \frac{m}{x}\rfloor\F(d)=\sum_{d|x}g(x)\mu(\frac{x}{d})=\sum_{d|x}\lfloor \frac{n}{x}\rfloor\lfloor \frac{m}{x}\rfloor\mu(\frac{x}{d})\\]
設\(T=\frac{x}{d}\)
\[ F(d)=\sum_{T=1}^{min(n,m)/d}\lfloor \frac{n}{Td}\rfloor\lfloor \frac{m}{Td}\rfloor\mu(T)\\]
那麽答案就是
\[ ans=\sum_{d=1}^{min(n,m)}d^{k}F(d)\ans=\sum_{d=1}^{min(n,m)}d^{k}\sum_{T=1}^{min(n,m)/d}\lfloor \frac{n}{Td}\rfloor\lfloor \frac{m}{Td}\rfloor\mu(T)\\]
然後這個式子已經可以預處理\(d^k\)，然後數論分塊來做了
不過會\(TLE\)，此時時間復雜度為\(O(nlogn+Tn)\)
所以我們還需要優化
考慮設\(k=Td\)
\[
ans=\sum_{k=1}^{min(n,m)}\lfloor \frac{n}{k}\rfloor\lfloor \frac{m}{k}\rfloor\sum_{d|k}\mu(\frac{k}{d})d^k
\]
考慮設\(s(d)=\sum_{d|k}\mu(\frac{k}{d})d^k\)
可以發現\(s\)是積性函數，因為\(\mu\)和\(d^k\)都是積性函數，積性函數和積性函數的狄利克雷卷積也是積性函數
所以\(s\)可以通過線篩求出
考慮如何線篩，對於\(a\)與\(b\)不互質的情況下，直接用\(s(ab)=s(a)*s(b)\)
當\(x\)為質數時，除了\(\mu(1)\)以及\(\mu(x)\)不為\(0\)，其余都為\(0\)，那麽\(s(x)=\mu(1)*x^k+\mu(x)*1^k\)
考慮\(a\)有\(x\)個質因子\(prime_i\)，那麽對於\(prime_i\)的貢獻\(s(prime_{i}^{x})=\mu(1)*prime_i^x+\mu(prime_i)*prime_i^{x-1}\)
那麽對於\(s(a*prime_i)\)，將相當於整體乘上一個\(prime_i^k\)
然後就可以篩了，我講的有點不清楚，可以看代碼理解一下
#include<cstdio>
#include<iostream>
#include<algorithm>
using namespace std;
void read(int &x) {
char ch; bool ok;
for(ok=0,ch=getchar(); !isdigit(ch); ch=getchar()) if(ch=='-') ok=1;
for(x=0; isdigit(ch); x=x*10+ch-'0',ch=getchar()); if(ok) x=-x;
}
#define rg register
const int maxn=5e6+10,mod=1e9+7;bool vis[maxn];
int T,k,n,m,pri[maxn],tot,ans,f[maxn],g[maxn];
int mi(int a,int b)
{
int ans=1;
while(b)
{
if(b&1)ans=1ll*ans*a%mod;
b>>=1,a=1ll*a*a%mod;
}
return ans;
}
void prepare()
{
f[1]=1;
for(rg int i=2;i<=5e6;i++)
{
if(!vis[i])pri[++tot]=i,g[i]=mi(i,k),f[i]=(g[i]-1+mod)%mod;
for(rg int j=1;pri[j]*i<=5e6&&j<=tot;j++)
{
vis[pri[j]*i]=1;
if(!(i%pri[j])){f[i*pri[j]]=1ll*f[i]*g[pri[j]]%mod;break;}
else f[i*pri[j]]=1ll*f[i]*f[pri[j]]%mod;
}
}
for(rg int i=1;i<=5e6;i++)f[i]=(f[i]+f[i-1])%mod;
}
int main()
{
read(T),read(k),prepare();
while(T--)
{
read(n),read(m);if(n>m)swap(n,m);ans=0;
for(rg int i=1,j;i<=n;i=j+1)
{
j=min(n/(n/i),m/(m/i));
ans=(ans+1ll*(f[j]-f[i-1]+mod)%mod*(n/i)%mod*(m/i)%mod)%mod;
}
printf("%d\n",ans);
}
}bzoj4407:於神之怒加強版