bzoj4407 于神之怒加强版

Description

---

给下N,M,K.求

$$\sum_{i=1}^{n}\sum_{j=1}^{m}{gcd(i,j)}^k \pmod{10^9+7}$$
1<=N,M,K<=5000000
1<=T<=2000

Solution

---

又补上一个大坑
前面的步骤其实同理就直接放最后结果好了

$$ans=\sum_{T=1}^{n}\lfloor\frac{n}{T}\rfloor \lfloor\frac{m}{T}\rfloor\sum_{d|T}d^k\mu(\frac{T}{d})$$
设 $g[i]=i^k$， $f[T]=\sum_{d|T}d^k\mu(\frac{T}{d})$，f是狄利克雷卷积，g和f都是积性函数。我们对于质数暴力快速幂，质数的次幂就递推，剩下的线性筛就能筛出来。注意prime的数组要开够不然会爆，大约开n/logn

讲一下怎么递推f。重点是筛出$f(p^a)$这样质数次幂的取值。显然有$f(p^a)=\sum_{d|p^a}d^k\mu(\frac{p^a}{d})=\sum_{i=0}^{a}p^{ki}\mu(p^{a-i})=\sum_{i=0}^{a}g(p^i)\mu(p^{a-i})$
观察一波发现当$a-i>1$时对结果无贡献，那么只需要考虑$i\geq a-1$的情况，原式等价于$f(p^a)=g(p^a)-g(p^{a-1})$，这样xjb搞搞就出来了

bzoj的权限好贵啊。。

Code

---

#include <stdio.h>
#include <string.h>
#include <algorithm>
#define rep(i,st,ed) for (int i=st;i<=ed;++i)

typedef long long LL;
const LL MOD=1000000007;
const int N=5000005;

int prime[N/10],mu[N+1],low[N+1];
bool not_prime[N+1];
LL f[N+1],g[N+1];

LL ksm(LL x,LL dep) {
    if (dep==0) return 1;
    if (dep==1) return x;
    LL tmp=ksm(x,dep/2);
    if (dep&1) return tmp*tmp%MOD*x%MOD;
    return tmp*tmp%MOD;
}

void pre_work(int n,LL k) {
    f[1]=mu[1]=g[1]=low[1]=1;
    rep(i,2,n) {
        if (!not_prime[i]) {
            prime[++prime[0]]=i;
            mu[i]=-1; low[i]=i;
            g[i]=ksm(i,k);
            f[i]=g[i]-1;
        }
        for (int j=1;j<=prime[0]&&i*prime[j]<=n;j++) {
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0) {
                low[i*prime[j]]=low[i]*prime[j];
                mu[i*prime[j]]=0;
                g[i*prime[j]]=g[i]*g[prime[j]]%MOD;
                f[i*prime[j]]=f[i/low[i]]*f[low[i]*prime[j]]%MOD;
                if (i==low[i]) {
                    f[i*prime[j]]=(g[i*prime[j]]-f[i]-g[i/prime[j]])%MOD;
                }
                break;
            }
            low[i*prime[j]]=prime[j];
            g[i*prime[j]]=g[i]*g[prime[j]]%MOD;
            f[i*prime[j]]=f[i]*f[prime[j]]%MOD;
            mu[i*prime[j]]=-mu[i];
        }
    }
    rep(i,1,n) f[i]=(f[i]+f[i-1])%MOD;
}

void solve(LL n,LL m) {
    if (m<n) std:: swap(n,m);
    LL ans=0;
    for (LL i=1,j;i<=n;i=j+1) {
        j=std:: min(n/(n/i),m/(m/i));
        ans=(ans+(n/i)*(m/i)%MOD*(f[j]-f[i-1])%MOD)%MOD;
    }
    if (ans<0) ans=(ans+MOD);
    printf("%lld\n", ans);
}

int main(void) {
    LL T,k; scanf("%lld%lld",&T,&k);
    pre_work(N,k);
    while (T--) {
        LL n,m; scanf("%lld%lld",&n,&m);
        solve(n,m);
    }
    return 0;
}