【BZOJ4407】于神之怒加强版

Description

给下N,M,K.求

Input

输入有多组数据，输入数据的第一行两个正整数T,K，代表有T组数据，K的意义如上所示，下面第二行到第T+1行，每行为两个正整数N,M，其意义如上式所示。
Output

如题
Sample Input

1 2

3 3
Sample Output

20

HINT

1<=N,M,K<=5000000,1<=T<=2000

Source

命题人：成都七中张耀楠，鸣谢excited上传。

300年没写过反演的题,今天再做一道内心几乎是崩溃的(摔
感谢鸟神提供的答疑w

题目要求求

$$\sum_{i=1}^n\sum_{j=1}^m gcd(i,j)^k$$
即为 $$\sum_{d=1}^{min(n,m)}d^k \sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=d]$$
后面那个东西很眼熟?似乎在梦中见过的样子w
$$\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=d]$$
$$=\sum_{i=1}^{\lfloor \frac n d\rfloor}\sum_{j=1}^{\lfloor \frac m d\rfloor}[gcd(i,j)=1]$$
显然可以反演咯w
此时的原式 $$\sum_{d=1}^{min(n,m)}d^k\sum_{i=1}^{\lfloor \frac n d\rfloor}\sum_{j=1}^{\lfloor \frac m d\rfloor}[gcd(i,j)=1]$$
$$=\sum_{d=1}^{min(n,m)}d^k\sum_{x=1}^{min(n,m)}\mu(x) \lfloor \frac n {xd}\rfloor \lfloor \frac m {xd}\rfloor$$
化到这里求解的复杂度是 $O(Tn)$,过不了.
怎么降复杂度? 丢掉d!
令 $y=xd$,现在d是y的约数了.取代一下.
原式
$$\sum_{d=1}^{min(n,m)}{d}^k\sum_{x=1}^{min(n,m)}\mu(\frac y d) \lfloor \frac n {y}\rfloor \lfloor \frac m {y}\rfloor$$
$$=\sum_{y=1}^{min(n,m)}\lfloor \frac n {y}\rfloor \lfloor \frac m {y}\rfloor \sum_{d|y}^{min(n,m)}\mu(\frac y d)d^k$$
里面那个是函数 $f(d)=d^k$与莫比乌斯函数的卷积.有前缀和就行了.
$d^k$这个可以放线筛里求,质数的话暴力快速幂,非质数的话就可以用两个因子的值搞出来.这样是 $O(n+T\sqrt n)$的.
如果直接对每个数暴力求得话是 $O(nlogn+T\sqrt n)$的.(跑50+s的那些都是这种)

最后我竟然跪在不会线性的求这个前缀和上了(摔

蓝月亮跑rank1好厉害呀

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define MAXN 5000010
#define GET (ch>='0'&&ch<='9')
#define P 1000000007
#define LL long long
using namespace std;
int T,k,n[2010],m[2010],maxn;
bool not_prime[MAXN];
int prime[MAXN],top;
int mu[MAXN],fac[MAXN];
LL f[MAXN],F[MAXN];
void in(int &x)
{
    char ch=getchar();x=0;
    while (!GET)    ch=getchar();
    while (GET) x=x*10+ch-'0',ch=getchar();
}
inline LL Pow(LL a,LL b)
{
    a%=P;LL ret=1;
    for (;b;b>>=1,a=a*a%P)  if (b&1)    ret*=a,ret%=P;
    return (ret+P)%P;
}
void check()
{
    mu[1]=1;F[1]=1;
    for (int i=2;i<=maxn;i++)
    {
        if (!not_prime[i])  prime[++top]=i,fac[i]=i,mu[i]=-1,f[i]=Pow(i,k),F[i]=f[i]-1;
        for (int j=1;j<=top&&i*prime[j]<=maxn;j++)
        {
            not_prime[i*prime[j]]=1;mu[i*prime[j]]=-mu[i];
            if (i%prime[j]==0)  
            {   
                mu[i*prime[j]]=0;fac[i*prime[j]]=fac[i]*prime[j];
                F[i*prime[j]]=fac[i]!=i?F[i/fac[i]]*F[fac[i]*prime[j]]%P:F[i]*f[prime[j]]%P;
                break;
            }
            F[i*prime[j]]=F[i]*F[prime[j]]%P;fac[i*prime[j]]=prime[j];
        }
    }
    for (int i=1;i<=maxn;i++)   mu[i]+=mu[i-1],F[i]+=F[i-1],F[i]%=P;
}
int main()
{
    int i;
    for (in(T),in(k),i=T;i;i--) in(n[i]),in(m[i]),maxn=max(maxn,max(n[i],m[i]));
    for (check();T;T--)
    {
        int N=n[T],M=m[T];
        LL ans=0;int t=min(N,M),last=1;
        for (int i=1;i<=t;i=last+1)
        {
            last=min(N/(N/i),M/(M/i));
            ans+=(F[last]-F[i-1]+P)*(N/i)%P*(M/i)%P;ans%=P;
        }
        printf("%lld\n",ans);
    }
}