【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∑n−1i=1∑n−1j=1(n−i)(n−j)[gcd(i,j)==1]+2n2−2n

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

ans′=∑n−1d=1μ[d]⌊(n−1)d⌋2n2−∑n−1d=1μ[d]dn⌊n−1d⌋2(1+⌊n−1d⌋)+∑n−1d=1μ[d]d2⌊n−1d⌋2(1+⌊n−1d⌋)24

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

对于∑n−1d=1μ[d]，s(n)=1−∑ni=2s(n/i)预处理一部分，哈希另一部分，使用分块可以快速求值。

同理对于∑n−1d=1μ[d]d，s(n)=1−∑ni=2is(n/i)

对于∑n−1d=1μ[d]d2，s(n)=1−∑ni=2i2s(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;
}