【jzoj5069】【GDSOI2017第二轮模拟】【蛋糕】【莫比乌斯反演】【杜教筛】

题目大意

CJY很喜欢吃蛋糕，于是YJC弄到了一块蛋糕，现在YJC决定和CJY分享蛋糕。

这块蛋糕上有n^2颗葡萄干，排成了一个n*n的点阵，每颗葡萄干互不相同且被编号为1~n^2。YJC决定沿着一条直线把蛋糕切成两份。YJC和CJY都很喜欢吃葡萄干，所以切出的两份蛋糕必须都包含至少一颗葡萄干。同时他们都不希望吃到不完整的葡萄干，所以切的时候不能经过任意一颗葡萄干。CJY喜欢1号葡萄干，所以他选择了包含1号葡萄干的那块蛋糕，而YJC则选择了另一块。

在吃蛋糕之前,CJY想知道有多少种切蛋糕的方案。两种方案是不同的当且仅当存在一颗葡萄干在一种方案中被分给了CJY，在另一种方案中被分给了YJC。CJY发现自己不会做这道题，所以他来向你请教。这个答案可能很大，你只需要告诉他答案mod p的结果就可以了。

解题思路

$ans=2\sum_{i=1}^{n-1}\sum_{j=1}^{n-1}(n-i)(n-j)[gcd(i,j)==1]+2n^2-2n$

对于前一部分通过莫比乌斯反演可得

$ans'=\sum_{d=1}^{n-1}\mu[d]\lfloor{(n-1)\over d}\rfloor^2n^2-\sum_{d=1}^{n-1}\mu[d]dn\lfloor{n-1\over d}\rfloor^2(1+\lfloor{n-1\over d}\rfloor)+\sum_{d=1}^{n-1}\mu[d]d^2{\lfloor{n-1\over d}\rfloor^2(1+\lfloor{n-1\over d}\rfloor)^2\over 4}$

我们对于带除法的部分可以分块处理，我们要求带$\mu$的前缀和，我们有杜教筛。

对于$\sum_{d=1}^{n-1}\mu[d]$，$s(n)=1-\sum_{i=2}^ns(n/i)$预处理一部分，哈希另一部分，使用分块可以快速求值。

同理对于$\sum_{d=1}^{n-1}\mu[d]d$，$s(n)=1-\sum_{i=2}^nis(n/i)$

对于$\sum_{d=1}^{n-1}\mu[d]d^2$，$s(n)=1-\sum_{i=2}^ni^2s(n/i)$

这样我们可以快速解决问题。

code

#include<set>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define LD double
#define LL long long
#define ULL unsigned long long
#define min(a,b) ((a<b)?a:b)
#define max(a,b) ((a>b)?a:b)
#define fo(i,j,k) for(int i=j;i<=k;i++)
#define fd(i,j,k) for(int i=j;i>=k;i--)
#define fr(i,j) for(int i=begin[j];i;i=next[i])
using namespace std;
int const mn=1e5+9,mm=2*1e5+9,lim=1e7,inf=1e9;
int n,p,ss[lim+9],tag[lim+9],mu[lim+9],mu2[lim+9],mu3[lim+9],hs[lim+9],a[lim+9],b[lim+9],c[lim+9];
int hash(int x){
    int tmp=x%lim;
    while(hs[tmp]&&(hs[tmp]!=x))tmp=(tmp+1==lim)?0:tmp+1;
    return tmp;
}
int s(int x){
    if(x<=lim)return mu[x];
    int ans=1,tmp=hash(x);
    if(hs[tmp]&&a[tmp])return a[tmp];
    for(int i=2,j;i<=x;i=j+1){
        j=x/(x/i);
        ans=(ans-1ll*s(x/i)*(j-i+1))%p;
    }
    hs[tmp]=x;
    return a[tmp]=ans;
}
int s2(int x){
    if(x<=lim)return mu2[x];
    int ans=1,tmp=hash(x);
    if(hs[tmp]&&b[tmp])return b[tmp];
    for(int i=2,j;i<=x;i=j+1){
        j=x/(x/i);
        ans=(ans-(1ll*(i+j)*(j-i+1)/2%p)*s2(x/i))%p;
    }
    hs[tmp]=x;
    return b[tmp]=ans;
}
int f2(int x){
    LL tmp=1ll*x*(x+1)/2;
    if(tmp%3==0)return tmp/3%p*(2*x+1)%p;
    else return 1ll*(2*x+1)/3*(tmp%p)%p;
}
int s3(int x){
    if(x<=lim)return mu3[x];
    int ans=1,tmp=hash(x);
    if(hs[tmp]&&c[tmp])return c[tmp];
    for(int i=2,j;i<=x;i=j+1){
        j=x/(x/i);
        ans=(ans-1ll*(f2(j)-f2(i-1))*s3(x/i))%p;
    }
    hs[tmp]=x;
    return c[tmp]=ans;
}
int f(int x){
    return (1ll*x*(x+1)/2%p)*(1ll*x*(x+1)/2%p)%p;
}
int main(){
    //freopen("cake.in","r",stdin);
    //freopen("cake.out","w",stdout);
    freopen("d.in","r",stdin);
    freopen("d.out","w",stdout);
    scanf("%d%d",&n,&p);
    mu[1]=1;
    mu2[1]=1;
    mu3[1]=1;
    fo(i,2,lim){
        if(!tag[i])ss[++ss[0]]=i,mu[i]=-1,mu2[i]=-i,mu3[i]=1ll*(-i)*i%p;
        fo(j,1,ss[0]){
            if(i*ss[j]>lim)break;
            tag[i*ss[j]]=1;
            mu[i*ss[j]]=-mu[i];
            mu2[i*ss[j]]=1ll*(-mu2[i])*ss[j]%p;
            mu3[i*ss[j]]=1ll*(-mu3[i])*ss[j]%p*ss[j]%p;
            if(i%ss[j]==0){
                mu[i*ss[j]]=0;
                mu2[i*ss[j]]=0;
                mu3[i*ss[j]]=0;
                break;
            }
        }
        mu[i]=(mu[i]+mu[i-1])%p;
        mu2[i]=(mu2[i]+mu2[i-1])%p;
        mu3[i]=(mu3[i]+mu3[i-1])%p;
    }
    int ans=0;
    for(int i=1,j;i<n;i=j+1){
        j=(n-1)/((n-1)/i);
        int tmp1=1ll*(s(j)-s(i-1))*((n-1)/i)%p*((n-1)/i)%p*n%p*n%p,
        tmp2=1ll*(s2(j)-s2(i-1))*n%p*((n-1)/i)%p*((n-1)/i)%p*((n-1)/i+1)%p,
        tmp3=1ll*(s3(j)-s3(i-1))*f((n-1)/i)%p%p;
        ans=(ans+1ll*(s(j)-s(i-1))*((n-1)/i)%p*((n-1)/i)%p*n%p*n%p
        -1ll*(s2(j)-s2(i-1))*n%p*((n-1)/i)%p*((n-1)/i)%p*((n-1)/i+1)%p
        +1ll*(s3(j)-s3(i-1))*f((n-1)/i)%p)%p;
    }
    printf("%d",((ans*2+1ll*2*n*n-2*n)%p+p)%p);
    return 0;
}