BZOJ 3309 莫比乌斯反演

题目链接 https://www.lydsy.com/JudgeOnline/problem.php?id=3309

题意：定义f(n)为n所含质因子的最大幂指数，求 $Ans=\sum _{i=1}^{a}\sum_{j=i}^{b}f(gcd(i,j))$

T<=10000

1<=a,b<=10^7

解析：考虑a<b

枚举最大公约数d，得到：

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

右边就是求$i\in[1,{\lfloor\frac{a}{d}\rfloor}],j\in[1,{\lfloor\frac{b}{d}\rfloor}]$ gcd(i,j)==1 的对数，根据莫比乌斯反演可以得到

$$Ans=\sum_{d=1}^a{f(d)}\sum_{d'=1}^{\lfloor\frac{a}{d}\rfloor}\mu(d'){\lfloor\frac{a}{dd'}\rfloor}{\lfloor\frac{b}{dd'}\rfloor}$$

令 T=dd'

$$Ans=\sum_{T=1}^{a}{\lfloor\frac{a}{T}\rfloor}{\lfloor\frac{b}{T}\rfloor}\sum_{d|T}\mu(\frac{T}{d}){f(d)}$$

左边可以分块求出来，所以右边要预处理出来前缀和，但是nlog(n) 的复杂度会超时，所以要大力分析一波

考虑 $G(n)=\sum_{d|n}\mu(\frac{n}{d}){f(d)}$ ，只有当 n/d是互不相等的素数乘积的形式才会产生贡献，其余都是0 ，

设 n=p1^q1*p2^q2*p3^q3...pi^qi , k=max{qi}, f(d) 的取值只有两种 k 或 k-1

考虑 n/d 的组合情况令，qi=k 的p集合为 A 剩余p集合为 B 当A集合全部都取时 f(d) = k-1 其余情况f(d) = k

1.当B不为空的时候 B集合的组合情况是奇数，偶数各占一半的。所以集合A 中的每一种情况与B组合也是奇偶数量相同，相互抵消了，G(n)=0。

2.当B为空的时候，也就是所有的qi是相同的都等于k，因为只有当A集合全部都取时 f(d) = k-1 ，假设A集合的大小为m，一共有2^m个组合方案：

1）m为奇数时 $\mu(\frac{n}{d})$ 是负的，f(d)= k 的组合方案数就是偶数多一个 G(n)=k-(k-1)=1

2) m为偶数时 $\mu(\frac{n}{d})$ 是正的，f(d)= k 的组合方案数就是奇数多一个 G(n)=(k-1)-k=-1

总结一下就是:

当n的素因子个数为偶数且素因子的幂次都相等时 G(n)=-1，

当n的素因子个数为奇数且素因子的幂次都相等时 G(n)=-1，

其余G(n)=0。

对于幂次等于1的情况我们可以线性筛求出来G(i)，G(i^k) =G(i) 指数级别的增长 O(n) 可以求出来所有的。

AC代码

#include <bits/stdc++.h>

#define pb push_back

#define mp make_pair

#define fi first

#define se second

#define all(a) (a).begin(), (a).end()

#define fillchar(a, x) memset(a, x, sizeof(a))

#define huan printf("\n");

#define debug(a,b) cout<<a<<" "<<b<<" ";

using namespace std;

const int maxn=1e7+,inf=0x3f3f3f3f;

typedef long long ll;

const ll mod = ;

typedef pair<int,int> pii;

int check[maxn],prime[maxn],G[maxn],sum[maxn];

void Mobius(int N)//线性筛求G(i)

{

    int pos=;G[]=;

    for (int i =  ; i <= N ; i++)

    {

        if (!check[i])

            prime[pos++] = i,G[i]=;

        for (int j =  ; j < pos && i*prime[j] <= N ; j++)

        {

            check[i*prime[j]] = ;

            if (i % prime[j] == )

            {

                G[i*prime[j]]=;

                break;

            }

            G[i*prime[j]]=-G[i];

        }

    }

}

int main()

{

    Mobius(1e7);

    sum[]=sum[]=;

    for(int i=;i<=1e7;i++)   //求G(i^k)

    {

        if(G[i]!=)

            for(ll j=i;j<=1e7;j*=i)

                sum[j]=G[i];

        sum[i]+=sum[i-];  //前缀和

    }

    int t,n,m;

    scanf("%d",&t);

    while(t--)

    {

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

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

        ll ans=;

        for(int i=,j;i<=n;i=j+) // 分块sqrt复杂度求出答案

        {

            j=min(n/(n/i),m/(m/i));

            ans+=1ll*(sum[j]-sum[i-])*(n/i)*(m/i);

        }

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

    }

}

巴特西

BZOJ 3309 莫比乌斯反演

最新文章

热门文章