bzoj 4407

莫比乌斯反演

还是推式子：

设$f(n)=n^{k}$

那就是上一道题了

推的过程如下：

$\sum_{i=1}^{a}\sum_{j=1}^{b}f(gcd(i,j))$

$\sum_{i=1}^{a}\sum_{j=1}^{b}\sum_{d=1}^{min(a,b)}[gcd(i,j)\equiv d]f(d)$

$\sum_{d=1}^{min(a,b)}f(d)\sum_{i=1}^{a}\sum_{j=1}^{b}[gcd(i,j)\equiv d]$

$\sum_{d=1}^{min(a,b)}f(d)\sum_{i=1}^{\frac{a}{d}}\sum_{j=1}^{\frac{b}{d}}[gcd(i,j)\equiv 1]$

$\sum_{d=1}^{min(a,b)}f(d)\sum_{i=1}^{\frac{a}{d}}\sum_{j=1}^{\frac{b}{d}}\sum_{t=1}^{min(\frac{a}{d},\frac{b}{d})}\mu(t)$

$\sum_{d=1}^{min(a,b)}f(d)\sum_{t=1}^{min(a,b)}\mu(t)\frac{a}{dt}\frac{b}{dt}$

令$T=dt$，得到：

$\sum_{T=1}^{min(a,b)}\frac{a}{T}\frac{b}{T}\sum_{d|T}f(d)\mu(\frac{T}{d})$

也就是：

$\sum_{T=1}^{min(a,b)}\frac{a}{T}\frac{b}{T}\sum_{d|T}d^{k}\mu(\frac{T}{d})$

考虑线性筛后面那堆东西，仍然分类讨论：

①.筛到的$p$与$i$互质：

此时我们考虑增加一个$p$的贡献，如果增加到$\mu$里，则原先那些直接取反

如果增加到$d^{k}$里，则相当于原先那些乘$p^{k}$

因此$g(ip)=(p^{k}-1)g(i)$

②.筛到的$p$与$i$不互质：

此时我们考虑增加一个$p$的贡献，如果增加到$\mu$里，则原先那些仍然取反

如果增加到$f$里，则原先那些多一个$p^{k}$的贡献

可...等等！

还有一种可能！

再考虑如果原先$\mu$里有一个$p$，然后增加到$f$里，此时会抵消掉取反的效果！

因此只需乘一个$p^{k}$即可

贴代码：

#include <cstdio>

#include <cmath>

#include <cstring>

#include <cstdlib>

#include <iostream>

#include <algorithm>

#include <queue>

#include <stack>

#define ll long long

using namespace std;

const ll mode=1000000007;

int mu[5000005];

int pri[5000005];

ll f[5000005];

ll pow_mul(ll x,ll y)

{

    ll ret=1;

    while(y)

    {

        if(y&1)ret=ret*x%mode;

        x=x*x%mode,y>>=1;

    }

    return ret;

}

bool used[10000005];

int cnt=0;

ll T,x,y,k;

void init()

{

    mu[1]=1;

    f[1]=1;

    for(int i=2;i<=5000000;i++)

    {

        if(!used[i])mu[i]=-1,pri[++cnt]=i,f[i]=(pow_mul(i,k)+mode-1)%mode;

        for(int j=1;j<=cnt&&i*pri[j]<=5000000;j++)

        {

            used[i*pri[j]]=1;

            if(i%pri[j]==0)

            {

                mu[i*pri[j]]=0;

                f[i*pri[j]]=f[i]*(f[pri[j]]+1)%mode;

                break;

            }

            mu[i*pri[j]]=-mu[i],f[i*pri[j]]=f[i]*f[pri[j]]%mode;

        }

    }

    for(int i=2;i<=5000000;i++)f[i]+=f[i-1],f[i]%=mode;

}

ll solve(ll a,ll b)

{

    ll las=1,ans=0;

    for(int i=1;i<=a&&i<=b;i=las+1)

    {

        las=min(a/(a/i),b/(b/i));

        ans+=(f[las]-f[i-1]+mode)*(a/i)%mode*(b/i)%mode;

        ans%=mode;

    }

    return ans;

}

template <typename T>inline void read(T &x)

{

    T f=1,c=0;char ch=getchar();

    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}

    while(ch>='0'&&ch<='9'){c=c*10+ch-'0';ch=getchar();}

    x=c*f;

}

int main()

{

    read(T),read(k);

    init();

    while(T--)

    {

        read(x),read(y);

        printf("%lld\n",solve(x,y));

    }

    return 0;

}

巴特西

bzoj 4407

最新文章

热门文章