洛谷

https://www.luogu.org/problemnew/show/P4449

(F(n)=sumlimits_{i=1}^{n}sumlimits_{i=1}^{m} gcd(i,j)^k)

首先加方括号,枚举g,提g:((min)表示(min(n,m))
(sumlimits_{g=1}^{min} g^k sumlimits_{i=1}^{n} sumlimits_{i=1}^{m} [gcd(i,j)==g])

(sumlimits_{g=1}^{min} g^k sumlimits_{i=1}^{lfloorfrac{n}{g} floor} sumlimits_{i=1}^{lfloorfrac{m}{g} floor} [gcd(i,j)==1])

后面莫比乌斯反演:(k被你用了真恶心)
(sumlimits_{g=1}^{min} g^k sumlimits_{x=1}^{min} mu(x){lfloorfrac{n}{gx} floor} {lfloorfrac{m}{gx} floor})

众所周知,这种情况要枚举(T=gx)
(sumlimits_{T=1}^{min}sumlimits_{g|T} g^k mu(frac{T}{g}) {lfloorfrac{n}{T} floor} {lfloorfrac{m}{T} floor})

提T:
(sumlimits_{T=1}^{min} {lfloorfrac{n}{T} floor} {lfloorfrac{m}{T} floor} sumlimits_{g|T} g^k mu(frac{T}{g}))

所以:
(F(n)=sumlimits_{i=1}^{n}sumlimits_{i=1}^{m} gcd(i,j)^k = sumlimits_{T=1}^{min} {lfloorfrac{n}{T} floor} {lfloorfrac{m}{T} floor} sumlimits_{g|T} g^k mu(frac{T}{g}))


(n,m)这么小那我给你搞个线性筛吧。记 (G(n)=sumlimits_{g|n} g^k mu(frac{n}{g})) ,这个可以 (O(nlogn)) 筛出来。但是我偏偏要线性筛。

每次回答一个分块了事了。


小心取模,靠。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

inline int read() {
    int x=0;
    char c=getchar();
    while(c<'0'||c>'9')
        c=getchar();
    do {
        x=(x<<3)+(x<<1)+c-'0';
        c=getchar();
    } while(c>='0'&&c<='9');
    return x;
}

inline void write(ll x) {
    if(x>9) {
        write(x/10);
    }
    putchar(x%10+'0');
    return;
}

const int MAXN=5e6;
const int MOD=1e9+7;

int pri[MAXN+1];
int &pritop=pri[0];
ll G[MAXN+1];
int pk[MAXN+1];

int k;

inline ll qpow(ll x,int n) {
    ll res=1;
    while(n) {
        if(n&1) {
            res*=x;
            if(res>=MOD)
                res%=MOD;
        }
        x*=x;
        if(x>=MOD)
            x%=MOD;
        n>>=1;
    }
    return res;
}

void sieve(int n=MAXN) {
    pk[1]=1;
    G[1]=1;
    for(int i=2; i<=n; i++) {
        if(!pri[i]) {
            pri[++pritop]=i;
            pk[i]=i;
            G[i]=qpow(i,k)-1ll;
            if(G[i]<0)
                G[i]+=MOD;
        }
        for(int j=1; j<=pritop; j++) {
            int &p=pri[j];
            int t=i*p;
            if(t>n)
                break;
            pri[t]=1;
            if(i%p) {
                pk[t]=pk[p];
                //积性函数
                G[t]=G[i]*G[p];
                if(G[t]>=MOD)
                    G[t]%=MOD;
            } else {
                pk[t]=pk[i]*p;
                if(pk[t]==t) {
                    //t是质数的幂次
                    G[t]=qpow(p,k)*G[i];
                    if(G[t]>=MOD)
                        G[t]%=MOD;
                } else {
                    //积性函数
                    G[t]=G[pk[t]]*G[t/pk[t]];
                    if(G[t]>=MOD)
                        G[t]%=MOD;
                }
                break;
            }
        }
    }

    for(int i=1;i<=n;i++){
        G[i]+=G[i-1];
        if(G[i]>=MOD)
            G[i]-=MOD;
    }
}

inline ll ans(int n,int m) {
    ll res=0;
    int N=min(n,m);
    for(int l=1,r;l<=N;l=r+1){
        r=min(n/(n/l),m/(m/l));
        ll tmp1=1ll*(n/l)*(m/l);
        if(tmp1>=MOD)
            tmp1%=MOD;
        ll tmp2=G[r]-G[l-1];
        if(tmp2<0)
            tmp2+=MOD;
        tmp1*=tmp2;
        if(tmp1>=MOD)
            tmp1%=MOD;
        res+=tmp1;
        if(res>=MOD)
            res%=MOD;
    }
    return res;
}

inline void solve() {
    int t=read();
    k=read();
    sieve();
    while(t--) {
        int n=read(),m=read();
        write(ans(n,m));
        putchar('
');
    }
}

int main() {
#ifdef Yinku
    freopen("Yinku.in","r",stdin);
#endif // Yinku
    solve();
    return 0;
}
原文地址:https://www.cnblogs.com/Yinku/p/10998146.html