BZOJ4407: 于神之怒加强版

BZOJ4407: 于神之怒加强版

Description

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

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


题解Here!

 

懵逼钨丝反演套路题。。。
题目要求:$$\sum_{i=1}^n\sum_{j=1}^mgcd(i,j)^k\mod (10^9+7)$$
先把那个$\mod (10^9+7)$丢一边去。。。
设$Ans=\sum_{i=1}^n\sum_{j=1}^mf(gcd(i,j)),f(x)=x^k$。
前面那个式子套路一下,可以看这题:BZOJ3529: [Sdoi2014]数表
于是:$$Ans=\sum_{D=1}^n\lfloor\frac{n}{D}\rfloor\lfloor\frac{m}{D}\rfloor\sum_{d|D}\mu(\frac{D}{d})f(d)$$
很显然嘛,前面的数论分块。
那$\sum_{d|D}\mu(\frac{D}{d})f(d)$呢?
你可能会说:暴力算到倍数里去!
$BUT$!这样是$O(n\log_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)\times f(b)$,所以$F(D)$也是积性函数。
质数时直接快速幂,质数的幂递推计算。
所以当$a,b$互质的时候,我们可以利用积性函数的性质$F(ab)=F(a)\times F(b)$。
关键是$a,b$不互质的时候怎么求解。
考虑:
$$F(p_i^{x_i})=\mu(1)\times (p_i^{x_i})^k+\mu(p_i)\times(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\n",solve(n,m));
    }
    return 0;
}

 

posted @ 2018-08-12 11:48  符拉迪沃斯托克  阅读(157)  评论(0编辑  收藏  举报
Live2D