【SDOI2014】数表【莫比乌斯反演】【树状数组】,电动牙刷产品设计

文件名：【SDOI2014】数表【莫比乌斯反演】【树状数组】,电动牙刷产品设计 【SDOI2014】数表【莫比乌斯反演】【树状数组】

传送门

传送门

题意：

$T$组询问给定$n,m,a$,求

$$\sum _{i=1}^N\sum _{j=1}^M[(\sum _{d\in N^{\ast }}[d\mid i][d\mid j]d)\le a]\sum _{d\in N^{\ast }}[d\mid i][d\mid j]d$$

$1\le n,m\le 1e5,1\le a\le 1e9,1\le T\le 2e4$

对$2^{31}$取模

反演神题

假设没有前面的限制

$$\sum _{i=1}^N\sum _{j=1}^M\sum _{d\in N^{\ast }}[d\mid i][d\mid j]$$

记$F(x)$表示约数和

$$\sum _{i=1}^N\sum _{j=1}^{M}F(gcd(i,j))$$

套路性地枚举$gcd$

$$\sum _d\sum _{i=1}^N\sum _{j=1}^M[gcd(i,j)=d]F(d)$$

套路性地提到前面

$$\sum _{d}F(d)\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=d]$$

后面是个套路性的反演

设

$$f(d)=\sum _{i=1}^N\sum _{j=1}^M[gcd(i,j)=d]$$

$$g(d)=\sum _{d\mid n}f(n)=\lfloor \frac{N}{d}\rfloor \lfloor \frac{M}{d}\rfloor$$

$$f(d)=\sum _{d\mid n}g(n)\mu (\frac{n}{d})=\sum _{d\mid n}\lfloor \frac{N}{n}\rfloor \lfloor \frac{M}{n}\rfloor \mu (\frac{n}{d})$$

代回去

$$\sum _{d}F(d)\sum _{d\mid n}\lfloor \frac{N}{n}\rfloor \lfloor \frac{M}{n}\rfloor \mu (\frac{n}{d})$$

$$\sum _{n=1}^{min(N,M)}\lfloor \frac{N}{n}\rfloor \lfloor \frac{M}{n}\rfloor \sum _{d\mid n}F(d)\mu (\frac{n}{d})$$

左边是个套路性的整除分块 不管

右边并没有很好的性质，只能筛出来暴力求……

等下，原题是有限制$F(d)\le a$的

设

$$G(n)=\sum _{d\mid n}F(d)\mu (\frac{n}{d})$$

我们发现只有不超过$a$的$F(d)$才对答案有贡献

而$T,N$都很小

所以可以把$F(d)$和$a$排个序，拿个指针指一下

跑到一个$a$把前面的$F(d)$枚举倍数统计更新$G$,再拿整除分块算答案

需要维护单点加，区间询问，用个树状数组即可

复杂度大概$O(NlogN+TlogT\sqrt{T})$

#include <iostream>#include <cstdio>#include <cstring>#include <cctype>#include <algorithm>#define MAXN 100005using namespace std;const int N=100000;int np[MAXN],pl[MAXN],cnt;int mu[MAXN],f[MAXN],mp[MAXN],pt[MAXN];inline int qpow(int a,int p){int ans=1;while (p){if (p&1) ans*=a;a*=a;p>>=1;}return ans;}void init(){np[1]=mu[1]=f[1]=1;for (int i=2;i<=N;++i){if (!np[i]){pl[++cnt]=i;mu[i]=-1;mp[i]=i;pt[i]=1;}int x;for (int j=1;(x=i*pl[j])<=N;++j){np[x]=1;mp[x]=pl[j];if (i%pl[j]==0){mu[x]=0;pt[x]=pt[i]+1;break;}mu[x]=-mu[i];pt[x]=1;}}for (int i=2;i<=N;i++){int t=qpow(mp[i],pt[i]);if (i==t) f[i]=f[i/mp[i]]+t;else f[i]=f[t]*f[i/t];}}int p[MAXN];inline bool cmp(const int& a,const int& b){return f[a]<f[b];}struct query{int n,m,a,id;}q[MAXN];inline bool operator <(const query& x,const query& y){return x.a<y.a;}struct BIT{int s[MAXN];inline int lowbit(const int& x){return x&-x;}inline void modify(int x,const int& v){for (;x<=N;s[x]+=v,x+=lowbit(x));}inline int query(int x){int ans=0;for (;x;ans+=s[x],x-=lowbit(x));return ans;}}g;int res[MAXN];//int calc(int x,int y,int a)//{//int ans=0;//for (int i=1;i<=x&&i<=y;i++)//if (x%i==0&&y%i==0)//ans+=i;//if (ans<=a) return ans;//else return 0;//}int main(){init();for (int i=1;i<=N;i++) p[i]=i;sort(p+1,p+N+1,cmp);int pos=0;int T;scanf("%d",&T);for (int i=1;i<=T;i++) {scanf("%d%d%d",&q[i].n,&q[i].m,&q[i].a);q[i].id=i;//int t=0;//for (int x=1;x<=q[i].n;x++)//for (int y=1;y<=q[i].m;y++)//t+=calc(x,y,q[i].a);//printf("%d\n",t&(~(1<<31)));}sort(q+1,q+T+1);for (int i=1;i<=T;i++){while (pos<N&&f[p[pos+1]]<=q[i].a){++pos;for (int d=1;d*p[pos]<=N;++d) g.modify(d*p[pos],f[p[pos]]*mu[d]);}if (q[i].n>q[i].m) swap(q[i].n,q[i].m);int ans=0;for (int l=1,r;l<=q[i].n;l=r+1){r=min(q[i].n/(q[i].n/l),q[i].m/(q[i].m/l));ans+=(q[i].n/l)*(q[i].m/l)*(g.query(r)-g.query(l-1));}res[q[i].id]=ans;}for (int i=1;i<=T;i++) printf("%d\n",res[i]&(~(1<<31)));return 0;}