P4449 于神之怒加强版

· · 题解

一个唯一一篇用杜教筛的题解。

题意：

求：

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

思路：

直接推式子：

\begin{aligned}\sum_{i=1}^n \sum_{j=1}^m \gcd(i,j)^k&=\sum_{t=1}^n \sum_{i=1}^n \sum_{j=1}^m [\gcd(i,j)=t] t^k\\&=\sum_{t=1}^n \sum_{i=1}^{\lfloor \frac{n}{t} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{t} \rfloor} [\gcd(i,j)=1] t^k \\&=\sum_{t=1}^n \sum_{i=1}^{\lfloor \frac{n}{t} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{t} \rfloor} \sum_{d|\gcd(i,j)} \mu(d) \times t^k \\ &= \sum_{t=1}^n \sum_{i=1}^{\lfloor \frac{n}{t} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{t} \rfloor} \sum_{d=1}^{\lfloor \frac{n}{t} \rfloor} [d|i][d|j]\mu(d) \times t^k \\&=\sum_{t=1}^n t^k \sum_{d=1}^{\lfloor \frac{n}{t} \rfloor} \mu(d) \sum_{i=1}^{\lfloor \frac{n}{t} \rfloor}[d|i] \sum_{j=1}^{\lfloor \frac{m}{t} \rfloor} [d|j] \\ &=\sum_{t=1}^n t^k \sum_{d=1}^{\lfloor \frac{n}{t} \rfloor} \mu(d) \times \lfloor \frac{n}{t \times d} \rfloor \times \lfloor \frac{m}{t \times d} \rfloor \end{aligned}

老套路，令 T = t \times d，原式化为：

\sum_{T=1}^n \lfloor \frac{n}{T} \rfloor \times \lfloor \frac{m}{T} \rfloor \sum_{t|T} t^k \times \mu(\frac{T}{k})

定义 k(x)=x^k，则：

\sum_{t|T} t^k \times \mu(\frac{T}{k})=(k \times \mu)(T)=\sum_{t|T} k(t) \times \mu(\frac{T}{k})

我们需要计算 (k \times \mu) 的前缀和，考虑杜教筛。

因为 \mu \times 1 = \varepsilon，则 k \times \mu \times 1=k \times \varepsilon，所以套入杜教筛板子：

s(n) = \sum_{i=1}^n(k \times \varepsilon)(i) - \sum_{i=2}^n s(\lfloor \frac{n}{i} \rfloor)

因为：

(k \times \varepsilon)(i)=\sum_{d|i}\varepsilon(d) \times k(\frac{i}{d})=k(i)

所以：

s(n)=\sum_{i=1}^n k(i) - \sum_{i=2}^n s(\lfloor \frac{n}{i} \rfloor)

预处理一下 k 的前缀和后即可计算。

杜教筛时只需要预处理 10^5 以内 (k \times \mu) 的前缀和即可，复杂度玄学，懒得证明了。

完整代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll N=5e6+10,K=1e5+10,mod=1e9+7;
inline ll read(){
    ll x=0,f=1;
    char c=getchar();
    while(c<'0'||c>'9'){
        if(c=='-')
          f=-1;
        c=getchar();
    }
    while(c>='0'&&c<='9'){
        x=(x<<1)+(x<<3)+(c^48);
        c=getchar();
    }
    return x*f;
}
inline void write(ll x){
    if(x<0){
        putchar('-');
        x=-x;
    }
    if(x>9)
      write(x/10);
    putchar(x%10+'0');
}
ll T,n,m,k,cnt=0; 
ll S[N],A[N],P[N],M[N];
bool f[N];
ll qpow(ll a,ll b){
    ll ans=1;
    while(b){
        if(b&1ll)
          ans=(ans*a)%mod;
        a=(a*a)%mod;
        b>>=1ll;
    }
    return ans;
}
void init(){
    M[1]=1;
    for(int i=2;i<N;i++){
        if(!f[i]){
            M[i]=-1;
            P[++cnt]=i;
        }
        for(ll j=1;j<=cnt&&i*P[j]<N;j++){
            f[i*P[j]]=1;
            if(i%P[j]==0)
              break;
            else 
              M[i*P[j]]=-M[i];
        }
    }
    for(int i=1;i<N;i++)
      S[i]=qpow(i,k);
    for(ll i=1;i<K;i++){
        if(!M[i])
          continue;
        for(ll j=i;j<K;j+=i)
          A[j]=(A[j]+M[i]*S[j/i]+mod)%mod;
    }
    for(int i=2;i<N;i++)
      S[i]=(S[i]+S[i-1])%mod;
    for(int i=2;i<K;i++)
      A[i]=(A[i]+A[i-1])%mod;
}
ll H(ll n){
    if(A[n])
      return A[n];
    ll ans=S[n],l,r;
    for(l=2;l<=n;l=r+1){
        r=n/(n/l);
        ans=(ans-((r-l+1)*H(n/l))%mod+mod)%mod;
    }
    A[n]=ans;
    return ans;
}
ll solve(ll n,ll m){
    ll ans=0,l,r;
    if(n>m)
      swap(n,m);
    for(l=1;l<=n;l=r+1){
        r=min(n/(n/l),m/(m/l));
        ans=(ans+((((n/l)*(m/l))%mod)*((H(r)-H(l-1)+mod)%mod))%mod)%mod;
    }
    return ans;
}
int main(){
    T=read(),k=read();
    init();
    while(T--){
        n=read(),m=read();
        write(solve(n,m));
        putchar('\n');
    }
    return 0;
}