P4449 于神之怒加强版

题目描述

给定 n,m,k ,计算

\sum_{i=1}^n \sum_{j=1}^m \gcd(i,j)^k

10^9 + 7 取模的结果。

Solution

Ans=\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k

Ans=\sum_{p=1}^np^k\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=p]

f(p)=\sum_{i=1}^n \sum_{j=1}^m [\gcd(i,j)=p]

F(p)=\sum_{p|d}f(d)=\sum_{p|d}\lfloor\frac{n}{d}\rfloor\times \lfloor \frac{m}{d} \rfloor

所以 f(p) = \sum_{p|n}\mu(\frac{n}{p})F(n)
Ans=\sum_{p=1}^np^k\sum_{p|d}\mu(\frac{d}{p})\lfloor\frac{n}{d}\rfloor\times \lfloor \frac{m}{d} \rfloor
更换枚举顺序
Ans=\sum_{d=1}^n\sum_{p|d}p^k\mu(\frac{d}{p})\lfloor\frac{n}{d}\rfloor\times \lfloor \frac{m}{d} \rfloor
g(x)=\sum_{p|d}=p^k\mu(\frac{d}{p})
这是一个狄利克雷卷积,也是积性函数,可以用线性筛筛出来

补充个线性筛处理积型函数的技巧

线性筛一般积性函数

若想线性筛出积性函数 f(x) ,就必须能够快速计算出以下函数值:

f(1)

f(p) p 是质数

f(p^k) p 是质数

其实就是含有的质因子个数小于等于 1 的所有数对应的函数值

常见的积性函数都会给出上述函数值的有关定义。对于自定义的一个积性函数(如狄利克雷卷积),就需要自行计算出上述函数值。

我们假设已经完成了上述函数值的计算,现在要求筛出所有至少含有两个质因子的数对应的函数值。
显然,一个含有至少两个质因子的数一定可以被分解成两个互质的且均不为 1 的数的乘积。此时我们就可以用 f(xy)=f(x)f(y) 计算得出相应的函数值。

我们考虑筛的过程中,i\times pri_j 会被 i 乘上 pri_j 筛掉

若将 i 唯一分解得到 p_1^{a_1}p_2^{a_2}\dots p_k^{a_k} ,那么一定有 pri_j\leq p_1

pri_jpri_ji 互质,可以直接 f(i\times pri_j)=f(i)\times f(pri_j)

pri_j=p_1 ,此时需要对 i 记录一个 low_i ,表示 i 中最小值因子的指数次幂,即 low_i=p_1^{a_1}

如果使用 i 除掉 low_i,那么结果中的最小值质因子一定大于 p_1 ,又因为 pri_j=p_1

从而可知 \gcd(\frac{i}{low_i},low_i\times pri_j)=1

于是 f(i\times pri_j)=f(\frac{i}{low_i})\times f(low_i\times pri_j)

void sieve()
{
    zhi[1]=low[1]=1;
    f[1]=对1直接定义;
    for (ll i=2;i<=n;i++)
    {
        if (!zhi[i]) low[i]=pri[++tot]=i,f[i]=对质数直接定义;
        for (ll j=1;j<=tot&&i*pri[j]<=n;j++)
        {
            zhi[i*pri[j]]=1;
            if (i%pri[j]==0)
            {
                low[i*pri[j]]=low[i]*pri[j];
                if (low[i]==i)
                    f[i*pri[j]]=对质数的若干次幂进行定义(一般由f[i]递推);
                else
                    f[i*pri[j]]=f[i/low[i]]*f[low[i]*pri[j]];
                break;
            }
            low[i*pri[j]]=pri[j];
            f[i*pri[j]]=f[i]*f[pri[j]];
        }
    }
}

于是这道题的代码:

#include<bits/stdc++.h>
#define ll long long 
const int N = 5e6+10;
const ll mod = 1e9+7;
int pri[N],vis[N],t,k,n,m;
ll low[N],g[N],cnt;
ll Mob[N],pre[N];
ll qpow(ll a,ll b){
    if(!b) return 1;
    ll c = qpow(a,b/2);
    c=c*c%mod;
    if(b&1) c=c*a%mod;
    return c;
}
void Mobi(){
    vis[1] = 1,g[1] = 1,low[1] = 1,Mob[1] = 1;
    for(int i=2;i<N;++i){
        if(!vis[i]){
            pri[++cnt] = i;
            Mob[i] = -1, low[i] = i , g[i] = (qpow(i,k)+Mob[i])%mod;
        }
        for(int j=1;j<=cnt&&1ll*pri[j]*i<N;++j){
            vis[i*pri[j]] = 1;
            if(i%pri[j]==0){
                low[i*pri[j]] = low[i]*pri[j];
                if(low[i]==i)
                    g[i*pri[j]]=(Mob[1]*qpow(i*pri[j],k)%mod+Mob[pri[j]]*qpow(i,k)%mod)%mod;
                else
                    g[i*pri[j]]=g[i/low[i]]*g[low[i]*pri[j]]%mod;
                break;
            }
            low[i*pri[j]] = pri[j];
            g[i*pri[j]] = g[i]*g[pri[j]]%mod;
        }
    }
    for(int i=1;i<N;++i)
        pre[i] = (pre[i-1]+g[i]) %mod;
}
int main(){
    scanf("%d%d",&t,&k);
    Mobi();
    while(t--){
        scanf("%d%d",&n,&m);
        if(n>m) std::swap(n,m);
        ll ans = 0 ;
        for(int l=1,r=0;l<=n;l=r+1){
            r = std::min(n,std::min(n/(n/l),m/(m/l)));
            ans = (ans+((pre[r]-pre[l-1])%mod)*(n/l)%mod*(m/l))%mod;
        }
        ans +=mod;
        ans%=mod;
        printf("%lld\n",ans);
    }
    return 0;
}

0 条评论

发表评论

Avatar placeholder

您的电子邮箱地址不会被公开。 必填项已用*标注