BZOJ4407: 于神之怒加强版

Description

给下N,M,K.求$$sum_{i=1}^nsum_{j=1}^mgcd(i,j)^kmod (10^9+7)$$

Input

输入有多组数据，输入数据的第一行两个正整数T,K，代表有T组数据，K的意义如上所示。

下面第二行到第T+1行，每行为两个正整数N,M，其意义如上式所示。

Output

如题

Sample Input

1 2
3 3

Sample Output

HINT

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

题解Here!

懵逼钨丝反演套路题。。。

题目要求：$$sum_{i=1}^nsum_{j=1}^mgcd(i,j)^kmod (10^9+7)$$

先把那个$mod (10^9+7)$丢一边去。。。

设$Ans=sum_{i=1}^nsum_{j=1}^mf(gcd(i,j)),f(x)=x^k$。

前面那个式子套路一下，可以看这题：BZOJ3529: [Sdoi2014]数表

于是：$$Ans=sum_{D=1}^nlfloorfrac{n}{D} floorlfloorfrac{m}{D} floorsum_{d|D}mu(frac{D}{d})f(d)$$

很显然嘛，前面的数论分块。

那$sum_{d|D}mu(frac{D}{d})f(d)$呢？

你可能会说：暴力算到倍数里去！

$BUT$！这样是$O(nlog_2n)$的，铁定$TLE$。。。

这次不一样了，$sum_{d|D}mu(frac{D}{d})f(d)$可以线性筛！

设$F(D)=sum_{d|D}mu(frac{D}{d})f(d)$。

我们发现函数$F(D)$是狄利克雷卷积的形式，并且$mu(i)$与$x^k$都是积性函数。

即：若$a,b$互质,$f(ab)=f(a) imes f(b)$，所以$F(D)$也是积性函数。
质数时直接快速幂，质数的幂递推计算。

所以当$a,b$互质的时候，我们可以利用积性函数的性质$F(ab)=F(a) imes F(b)$。

关键是$a,b$不互质的时候怎么求解。

考虑：

$$F(p_i^{x_i})=mu(1) imes (p_i^{x_i})^k+mu(p_i) imes(p_i^{x_{i-1}})^k$$

剩下的项$mu$都是$0$。

那么对于一个质因数来说，$F(p_i^{x_i})$变成$F(p_i^{x_{i+1}})$的影响就是在$F(p_i^{x_i})$的基础上乘$p_i^k$。

其它质因数与其互质，直接乘起来就好啦。

那么对于函数$F(D)$我们就可用线性筛进行预处理，搞个前缀和就好了。

具体的看代码。

附代码：

#include<iostream>
#include<algorithm>
#include<cstdio>
#define MAXN 5000010
#define MOD 1000000007LL
using namespace std;
int k=0,prime[MAXN];
long long p,f[MAXN],g[MAXN],sum[MAXN];
bool np[MAXN];
inline int read(){
    int date=0,w=1;char c=0;
    while(c<'0'||c>'9'){if(c=='-')w=-1;c=getchar();}
    while(c>='0'&&c<='9'){date=date*10+c-'0';c=getchar();}
    return date*w;
}
long long mexp(long long a,long long b,long long c){
    long long s=1;
    while(b){
        if(b&1)s=s*a%c;
        a=a*a%c;
        b>>=1;
    }
    return s;
}
void make(){
    int m=MAXN-10;
    f[1]=g[1]=1;
    for(int i=2;i<=m;i++){
        if(!np[i]){
            prime[++k]=i;
            g[i]=mexp(i,p,MOD);
            f[i]=g[i]-1;
        }
        for(int j=1;j<=k&&prime[j]*i<=m;j++){
            np[prime[j]*i]=true;
            if(i%prime[j]==0){
                f[prime[j]*i]=f[i]*g[prime[j]]%MOD;
                break;
            }
            f[prime[j]*i]=f[i]*f[prime[j]]%MOD;
        }
    }
    for(int i=1;i<=m;i++)sum[i]=(sum[i-1]+f[i])%MOD;
}
long long solve(long long n,long long m){
    long long ans=0;
    if(n>m)swap(n,m);
    for(int i=1,last=1;i<=n;i=last+1){
        last=min(n/(n/i),m/(m/i));
        ans=(ans+(sum[last]-sum[i-1]+MOD)%MOD*(n/i)%MOD*(m/i)%MOD)%MOD;
    }
    return ans;
}
int main(){
    int t=read();
    p=read();
    make();
    while(t--){
        int n=read(),m=read();
        printf("%lld
",solve(n,m));
    }
    return 0;
}