bzoj4407 于神之怒加强版

Description


给下N,M,K.求

i=1nj=1mgcd(i,j)k(mod109+7) ∑ i = 1 n ∑ j = 1 m g c d ( i , j ) k ( mod 10 9 + 7 )

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

Solution


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

ans=T=1nnTmTd|Tdkμ(Td) a n s = ∑ T = 1 n ⌊ n T ⌋ ⌊ m T ⌋ ∑ d | T d k μ ( T d )

g[i]=ik g [ i ] = i k f[T]=d|Tdkμ(Td) f [ T ] = ∑ d | T d k μ ( T d ) ,f是狄利克雷卷积,g和f都是积性函数。我们对于质数暴力快速幂,质数的次幂就递推,剩下的线性筛就能筛出来。注意prime的数组要开够不然会爆,大约开n/logn

讲一下怎么递推f。重点是筛出 f(pa) f ( p a ) 这样质数次幂的取值。显然有 f(pa)=d|padkμ(pad)=ai=0pkiμ(pai)=ai=0g(pi)μ(pai) f ( p a ) = ∑ d | p a d k μ ( p a d ) = ∑ i = 0 a p k i μ ( p a − i ) = ∑ i = 0 a g ( p i ) μ ( p a − i )
观察一波发现当 ai>1 a − i > 1 时对结果无贡献,那么只需要考虑 ia1 i ≥ a − 1 的情况,原式等价于 f(pa)=g(pa)g(pa1) 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;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值