bzoj 3561

题意：求$\sum_{i=1}^{n}\sum_{j=1}^{n}lcm(i,j)^{gcd(i,j)}$

神仙题...

首先可能会想到一个转化，就是$lcm(i,j)=\frac{ij}{gcd(i,j)}$

然后大力往下推式子，发现你推不下去了...

因为$d$在分母上！！！

然后我们考虑换一种推法：如果我们对$ij$同时除掉$gcd(i,j)$，这样的话问题就可以转化成这个样子：

$\sum_{d=1}^{n}\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}[gcd(i,j)\equiv 1](ijd)^{d}$

然后把$d$拿出来，维护一下后面那坨，有：

$\sum_{d=1}^{n}d^{d}\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}\sum_{t|gcd(i,j)}\mu(t)(ij)^{d}$

改变一下枚举顺序，得到：

$\sum_{d=1}^{n}d^{d}\sum_{t=1}^{\frac{n}{d}}\mu(t)\sum_{i=1}^{\frac{n}{dt}}\sum_{j=1}^{\frac{m}{dt}}(ijt^{2})^{d}$

（也就是在后面的$ij$乘积这一项中单独考虑$t$的贡献）

然后再整理，就得到：

$\sum_{d=1}^{n}d^{d}\sum_{t=1}^{\frac{n}{d}}\mu(t)t^{2d}\sum_{i=1}^{\frac{n}{dt}}i^{d}\sum_{j=1}^{\frac{m}{dt}}j^{d}$

然后我们就暴力计算即可，每次计算时都要先预处理出$[1,\frac{n}{d}]$的$i^{d}$的前缀和，再暴力查询即可，时间复杂度为调和级数$O(nlnn)$

注意每次求幂可以递推，不要快速幂！！！会退化成$O(nlog_{2}^{2}n)$！

代码：

#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 pri[5000005],mu[5000005],used[5000005];

ll S[5000005],mi[5000005];

int cnt=0;

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;

}

void init()

{

    mu[1]=1;

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

    {

        if(!used[i])pri[++cnt]=i,mu[i]=-1;

        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;break;}

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

        }

    }

}

int main()

{

    init();

    ll n,m;

    scanf("%lld%lld",&n,&m);

    if(n>m)swap(n,m);

    ll ans=0;

    for(int i=1;i<=m;i++)mi[i]=1;

    for(int i=1;i<=n;i++)

    {

        ll s=pow_mul(i,i);

        S[0]=0;

        for(int j=1;j<=m/i;j++)mi[j]=mi[j]*j%mode,S[j]=(S[j-1]+mi[j])%mode;

        ll tempc=0;

        for(int j=1;j<=n/i;j++)

        {

            ll temps=(mu[j]*mi[j]*mi[j]%mode+mode)%mode;

            tempc=(tempc+temps*S[n/i/j]%mode*S[m/i/j]%mode)%mode;

        }

        ans=(ans+s*tempc)%mode;

    }

    printf("%lld\n",ans);

    return 0;

}

巴特西

bzoj 3561

最新文章

热门文章