计蒜客 ICPC沈阳网络赛 Convex Hull（容斥原理 + 莫比乌斯函数）

本文链接：https://blog.csdn.net/u013534123/article/details/82695640

大致提议：求题目中所给式子的求和。

真的是非常非常非常简单的一道题目，式子都推对了，我竟然去想杜教筛而没有想容斥……

根据函数的定义，很容易发现这个函数其实是 $gay(i)=|\mu(i)|*i^2$ 或者说 $gay(i)=\mu(i)^2*i^2$ 。然后对于题目给的式子，我们可以交换一下求和次序，可以得到：

$\large ans=\sum_{i=1}^{n}(n-i+1)*gay(i)$

然后把它拆开成为两部分：

$\large \begin{align*} ans&=\sum_{i=1}^{n}(n+1)*gay(i)-\sum_{i=1}^{n}i*gay(i) \\ &=\sum_{i=1}^{n}\mu(i)^2*i^2-\sum_{i=1}^{n}\mu(i)^2*i^3\end{align*}$

比赛的时候推到这一步，然后就想着1e10的数据，肯定杜教筛之类的，然后就自闭了……

其实这里直接上容斥原理就行啦。仔细观察两个求和式子，当莫比乌斯函数值不为0的时候，就把i^2或者i^3加上。那么我在加的时候，可以考虑直接把它的所有倍数的加上，然后在容斥掉重复的。具体来说，当i为1的时候，把所有数字的平方和立方和加上，然后走到任意一个莫比乌斯函数值不为0的i的时候，把以及对应的立方都筛掉。你会发现，这个其实是：

$\large \begin{align*}\sum_{j=1}^{\lfloor\frac{x}{i^2}\rfloor} (j*i^2)^2&=i^4*\sum_{i=1}^{\lfloor\frac{x}{i^2}\rfloor}j^2 \\ &=i^4*\frac{\lfloor\frac{x}{i^2}\rfloor(\lfloor\frac{x}{i^2}\rfloor+1)*(2\lfloor\frac{x}{i^2}\rfloor+1)}{6}\end{align*}$

同理，对于另一个i^3的求和式子也是一样可以做同样的变换。所以说，我们可以快速的求出这两个部分的和，因此容斥原理直接做即可。具体见代码：

#include<bits/stdc++.h>
//#define mod 1000000007
#define LL long long
#define pb push_back
#define lb lower_bound
#define ub upper_bound
#define INF 0x3f3f3f3f
#define sf(x) scanf("%lld",&x)
#define sc(x,y,z) scanf("%d%d%d",&x,&y,&z)
#define clr(x,n) memset(x,0,sizeof(x[0])*(n+5))
#define file(x) freopen(#x".in","r",stdin),freopen(#x".out","w",stdout)
using namespace std;

const int N = 1e5 + 10;
const int len = 20;

int mu[N],p[N],sz;
bool isp[N];
LL n,mod;

inline LL multi(LL A,LL B,LL mod)
{
    LL a=A>>len,b=A&((1<<len)-1);
    LL c=B>>len,d=B&((1<<len)-1);
    return ((((((a*c)<<len)+b*c+a*d)%mod)<<len)+b*d)%mod;
}

void init()
{
    sz=0;
    mu[1]=1;
    for(int i=2;i<N;i++)
    {
        if(!isp[i])p[++sz]=i,mu[i]=-1;
         for(int j=1;j<=sz&&(LL)i*p[j]<N;j++)
         {
            isp[i*p[j]]=1;
            if(i%p[j]==0)
            {
                mu[p[j]*i]=0; break;
            } else mu[p[j]*i]=-mu[i];
        }
    }
}

inline LL s2(LL x)
{
    return multi(multi(x,x+1,mod*6),2*x+1,mod*6)/6;
}

inline LL s3(LL x)
{
    LL res=multi(x,x+1,mod<<1)>>1;
    return multi(res,res,mod);
}

int main()
{
    init();
    while(~sf(n))
    {
        sf(mod); LL res=0;
        for(int i=1;(LL)i*i<=n;i++)
        {
            if (!mu[i]) continue;
            LL t1=(LL)i*i;
            LL t2=multi(t1,t1,mod);
            res+=mu[i]*multi(s2(n/t1),multi(n+1,t2,mod),mod);
            res-=mu[i]*multi(s3(n/t1),multi(t2,t1,mod),mod);
            res%=mod;
        }
        printf("%lld\n",(res+mod)%mod);
    }
    return 0;
}