并不对劲的p4449于神之怒加强版

题目大意

给定(t,k(tleq2000,kleq5*10^6))
(t)组询问，每组给出(n,m(n,mleq5*10^6))求$sum_{i=1}^n sum_{j=1}^m mathrm{gcd}(i,j)^k $

题解

假设(n)较小
枚举gcd
原式=(sum_{a=1}^{n}a^ksum_{i=1}^{lfloorfrac{n}{a} floor} sum_{j=1}^{lfloorfrac{m}{a} floor} [gcd(i,j)=1])
=(sum_{a=1}^{n}a^ksum_{i=1}^{lfloorfrac{n}{a} floor} sum_{j=1}^{lfloorfrac{m}{a} floor}sum_{d|i且d|j}mu(d))
=(sum_{a=1}^{n}a^ksum_{d=1}^{lfloorfrac{n}{a} floor}mu(d)*{lfloorfrac{n}{d*a} floor}*{lfloorfrac{m}{d*a} floor})
设(b=d*a)
=(sum_{a=1}^{n}a^ksum_{a|b}^{n}mu(frac{b}{a})*{lfloorfrac{n}{b} floor}*{lfloorfrac{m}{b} floor})
=(sum_{b=1}^{n}{lfloorfrac{n}{b} floor}*{lfloorfrac{m}{b} floor}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('
');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('
');
	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;
}