洛谷 - P4449 - 于神之怒加强版

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

\(F(n)=\sum\limits_{i=1}^{n}\sum\limits_{i=1}^{m} gcd(i,j)^k\)

首先加方括号，枚举g，提g：（\(min\)表示\(min(n,m)\)）

\(\sum\limits_{g=1}^{min} g^k \sum\limits_{i=1}^{n} \sum\limits_{i=1}^{m} [gcd(i,j)==g]\)

\(\sum\limits_{g=1}^{min} g^k \sum\limits_{i=1}^{\lfloor\frac{n}{g}\rfloor} \sum\limits_{i=1}^{\lfloor\frac{m}{g}\rfloor} [gcd(i,j)==1]\)

后面莫比乌斯反演：（k被你用了真恶心）

\(\sum\limits_{g=1}^{min} g^k \sum\limits_{x=1}^{min} \mu(x){\lfloor\frac{n}{gx}\rfloor} {\lfloor\frac{m}{gx}\rfloor}\)

众所周知，这种情况要枚举\(T=gx\)：

\(\sum\limits_{T=1}^{min}\sum\limits_{g|T} g^k \mu(\frac{T}{g}) {\lfloor\frac{n}{T}\rfloor} {\lfloor\frac{m}{T}\rfloor}\)

提T：

\(\sum\limits_{T=1}^{min} {\lfloor\frac{n}{T}\rfloor} {\lfloor\frac{m}{T}\rfloor} \sum\limits_{g|T} g^k \mu(\frac{T}{g})\)

所以：

\(F(n)=\sum\limits_{i=1}^{n}\sum\limits_{i=1}^{m} gcd(i,j)^k = \sum\limits_{T=1}^{min} {\lfloor\frac{n}{T}\rfloor} {\lfloor\frac{m}{T}\rfloor} \sum\limits_{g|T} g^k \mu(\frac{T}{g})\)

\(n,m\)这么小那我给你搞个线性筛吧。记 \(G(n)=\sum\limits_{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('\n');

    }

}

int main() {

#ifdef Yinku

    freopen("Yinku.in","r",stdin);

#endif // Yinku

    solve();

    return 0;

}

巴特西

洛谷 - P4449 - 于神之怒加强版 - 莫比乌斯反演

最新文章

热门文章