并不对劲的p4449于神之怒加强版
题目大意
给定\(t,k(t\leq2000,k\leq5*10^6)\)
\(t\)组询问,每组给出\(n,m(n,m\leq5*10^6)\)求$\sum_{i=1}^n \sum_{j=1}^m \mathrm{gcd}(i,j)^k $
题解
假设\(n\)较小
枚举gcd
原式=\(\sum_{a=1}^{n}a^k\sum_{i=1}^{\lfloor\frac{n}{a}\rfloor} \sum_{j=1}^{\lfloor\frac{m}{a}\rfloor} [gcd(i,j)=1]\)
=\(\sum_{a=1}^{n}a^k\sum_{i=1}^{\lfloor\frac{n}{a}\rfloor} \sum_{j=1}^{\lfloor\frac{m}{a}\rfloor}\sum_{d|i且d|j}\mu(d)\)
=\(\sum_{a=1}^{n}a^k\sum_{d=1}^{\lfloor\frac{n}{a}\rfloor}\mu(d)*{\lfloor\frac{n}{d*a}\rfloor}*{\lfloor\frac{m}{d*a}\rfloor}\)
设\(b=d*a\)
=\(\sum_{a=1}^{n}a^k\sum_{a|b}^{n}\mu(\frac{b}{a})*{\lfloor\frac{n}{b}\rfloor}*{\lfloor\frac{m}{b}\rfloor}\)
=\(\sum_{b=1}^{n}{\lfloor\frac{n}{b}\rfloor}*{\lfloor\frac{m}{b}\rfloor}\sum_{a|b}^{n}\mu(\frac{b}{a})*a^k\)
设\(f(x)=\sum_{a|x}^{n}\mu(\frac{x}{a})*a^k\)
发现\(f(x)\)是\(\mu\)卷\(x^k\),所以\(f(x)\)也是个积性函数,可以用筛法求
要考虑两件事:1.当\(p\)为质数时,\(f(p)=?\) 2.当\(p\)为质数且\(p|q\)时,\(f(p*q)=?\)
1.\(f(p)=1^k*\mu(p)+p^k*mu(1)=p^k-1\)
2.设\(q\)中\(p\)这个因子的指数为\(c\),则有\(f(q)=f(\frac{q}{p^c})*f(p^c),f(p*q)=f(\frac{q}{p^c})*f(p^{c+1})\),即\(f(p*q)=f(q)*(\frac{f(p^{c+1})}{f(p^c)})\)
那么只要知道\(\frac{f(p^{c+1})}{f(p^c)}\),就能用\(f(q)\)推出\(f(p*q)\)了
发现\(f(p^c)=\sum_{a|p^c}a^k*\mu(\frac{p^c}{a})\),其中\(\mu(\frac{p^c}{a})\)只在\(a=p^c\)时为1,在\(a=p^{c-1}\)时为-1,其余时刻都为0
那就有\(f(p^c)=p^{c*k}-p^{(c-1)*k}\)
同理可得\(f(p^{c+1})=p^{(c+1)*k}-p^{c*k}\)
所以\(\frac{f(p^{c+1})}{f(p^c)}=p^k\)
预处理筛出\(f(x)\)后,询问时整出分块就行了
代码
#include<algorithm>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<ctime>
#include<iomanip>
#include<iostream>
#include<map>
#include<queue>
#include<set>
#include<stack>
#include<vector>
#define rep(i,x,y) for(register int i=(x);i<=(y);++i)
#define dwn(i,x,y) for(register int i=(x);i>=(y);--i)
#define maxn 5000010
#define lim (maxn-10)
#define LL long long
#define add(x) (x>=mod?x-mod:x)
using namespace std;
int read()
{
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)&&ch!='-')ch=getchar();
if(ch=='-')f=-1,ch=getchar();
while(isdigit(ch))x=(x<<1)+(x<<3)+ch-'0',ch=getchar();
return x*f;
}
void write(int x)
{
if(x==0){putchar('0'),putchar('\n');return;}
int f=0;char ch[20];
if(x<0)putchar('-'),x=-x;
while(x)ch[++f]=x%10+'0',x/=10;
while(f)putchar(ch[f--]);
putchar('\n');
return;
}
const LL mod=1000000007;
int no[maxn],p[maxn],cnt,kth[maxn],mu[maxn],f[maxn],n,m,k,t;
int mul(int x,int y){int res=1;while(y){if(y&1)res=(LL)res*(LL)x%mod;x=(LL)x*(LL)x%mod,y>>=1;}return res;}
int main()
{
t=read(),k=read();
no[1]=f[1]=1;
rep(i,2,lim)
{
if(!no[i])p[++cnt]=i,kth[i]=mul(i,k),f[i]=(-1+kth[i]+mod)%mod;
for(int j=1;j<=cnt&&i*p[j]<=lim;j++)
{
no[i*p[j]]=1;
if(i%p[j]==0){f[i*p[j]]=(LL)f[i]*(LL)kth[p[j]]%mod;break;}
f[i*p[j]]=(LL)f[p[j]]*(LL)f[i]%mod;
}
}
rep(i,1,lim)f[i]=add(f[i]+f[i-1]);
while(t--)
{
n=read(),m=read();if(n>m)swap(n,m);
int ans=0;
for(int l=1,r;l<=n;l=r+1)
{
r=min(n/(n/l),m/(m/l));
ans=(ans+(LL)(n/l)*(LL)(m/l)%mod*(LL)((f[r]-f[l-1]+mod)%mod)%mod)%mod;
}
write(ans);
}
return 0;
}