您的位置:首页 > 其它

3529: [Sdoi2014]数表|莫比乌斯反演|树状数组

2016-03-17 20:52 549 查看
以下式子均设n≤mn\leq m

先设M(x)=∑d|xdM(x)=\sum_{d|x}d

在不考虑a的限制下

Ans=∑i=1n∑j=1nM(gcd(i,j))Ans=\sum_{i=1}^{n}\sum_{j=1}^nM(gcd(i,j))

=∑i=1nM(i)∗sum(i)=\sum_{i=1}^{n}M(i)*sum(i)

sum(i)sum(i)表示gcd(x,y)=igcd(x,y)=i的(x,y)(x,y)的个数

显然这个反演一下就可以得到sum(i)=∑x=1⌊ni⌋u(i)∗⌊ni∗x⌋∗⌊mi∗x⌋sum(i)=\sum_{x=1}^{\lfloor\frac{n}{i}\rfloor}u(i)*\lfloor\frac{n}{i*x}\rfloor*\lfloor\frac{m}{i*x}\rfloor

Ans=∑i=1nM(i)∗∑x=1⌊ni⌋u(i)∗⌊ni∗x⌋∗⌊mi∗x⌋Ans=\sum_{i=1}^{n}M(i)*\sum_{x=1}^{\lfloor\frac{n}{i}\rfloor}u(i)*\lfloor\frac{n}{i*x}\rfloor*\lfloor\frac{m}{i*x}\rfloor

=∑T=1n⌊nT⌋∗⌊mT⌋∗(∑d|Tu(d)∗M(Td))=\sum_{T=1}^{n}\lfloor\frac{n}{T}\rfloor*\lfloor\frac{m}{T}\rfloor*\big(\sum_{d|T}u(d)*M(\frac{T}{d})\big)

这样只需要线性筛搞出后面的这个∑d|Tu(d)∗M(Td)\sum_{d|T}u(d)*M(\frac{T}{d})就可以n−√\sqrt n处理每次询问了

然后再回来考虑a的限制,只有M(x)≤aM(x)\leq a才会对答案有贡献,所以可以对询问按aa的大小排序,每次都把小于当前的a的M(x)M(x)加入到树状数组中,这样离线处理每次询问

然后处理模数时中途直接全部intint最后在andand一下intint的最大值就可以

#include<algorithm>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<vector>
#include<cstdio>
#include<queue>
#include<cmath>
#include<set>
#include<map>
#define ll long long
#define N 100005
using namespace std;
int sc()
{
int i=0,f=1; char c=getchar();
while(c>'9'||c<'0'){if(c=='-')f=-1;c=getchar();}
while(c>='0'&&c<='9')i=i*10+c-'0',c=getchar();
return i*f;
}
struct W{int n,m,p,a;}a
;
int low
,A
,prime
,tr
;
int M
,u
,ans
,pos
;
void pre()
{
int top=0;M[1]=u[1]=1;
for(int i=2;i<N;i++)
{
if(!A[i])
{
u[low[i]=prime[++top]=i]=-1;
M[i]=i+1;
}
for(int j=1;i*prime[j]<N;j++)
{
A[prime[j]*i]=1;
if(i%prime[j]==0)
{
low[i*prime[j]]=low[i]*prime[j];
if(low[i]==i)
M[i*prime[j]]=M[i]+i*prime[j];
else M[i*prime[j]]=M[i/low[i]]*M[low[i]*prime[j]];
break;
}
u[prime[j]*i]=-u[i];
low[prime[j]*i]=prime[j];
M[prime[j]*i]=M[i]*(prime[j]+1);
}
}
}
bool cmp(W a,W b){return a.a<b.a;}
bool Cmp(int a,int b){return M[a]<M[b];}
void change(int x,int f)
{
for(;x<N;x+=x&-x)tr[x]+=f;
}
int ask(int x)
{
int ans=0;
for(;x;x-=x&-x)ans+=tr[x];
return ans;
}
int query(int n,int m)
{
if(n>m)swap(n,m);
int ans=0;
for(int i=1,last;i<=n;i=last+1)
{
last=min(n/(n/i),m/(m/i));
ans+=(n/i)*(m/i)*(ask(last)-ask(i-1));
}
return ans;
}
int main()
{
pre();
int W=sc();
for(int i=1;i<=W;i++)
a[i].n=sc(),a[i].m=sc(),a[i].a=sc(),a[i].p=i;
sort(a+1,a+W+1,cmp);
for(int i=1;i<N;i++)pos[i]=i;
sort(pos+1,pos+N,Cmp);
int now=1;
for(int i=1;i<=W;i++)
{
while(now<N&&M[pos[now]]<=a[i].a)
{
int x=pos[now++];
for(int j=x;j<N;j+=x)
change(j,M[x]*u[j/x]);
}
ans[a[i].p]=query(a[i].m,a[i].n);
}
for(int i=1;i<=W;i++)
printf("%d\n",ans[i]&0x7FFFFFFF);
return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: