BZOJ3529 || 洛谷P3312 [SDOI2014]数表【莫比乌斯反演+树状数组】

Time Limit: 10 Sec
Memory Limit: 512 MB

Description

有一张 n×m 的数表，其第 i 行第 j 列（1 <= i <= n, 1 <= j <= m）的数值为
能同时整除 i 和 j 的所有自然数之和。给定 a , 计算数表中不大于 a 的数之和。

Input

输入包含多组数据。
输入的第一行一个整数Q表示测试点内的数据组数
接下来Q行，每行三个整数n，m，a(|a| < =10^9)描述一组数据。
1 < =N．m < =10^5 ， 1 < =Q < =2×10^4

Output

对每组数据，输出一行一个整数，表示答案模2^31的值

---

题目分析

若不考虑a，那么可以列出答案的表达式

$$ans=\sum _{i=1}^n\sum _{j=1}^n\sigma (gcd(i,j))$$

(其中$\sigma (k)$为$k$的约数和)

既然和gcd有关，应该马上想到有关gcd的莫比乌斯反演的常见套路
$g(k)$为给定范围内满足$gcd(x,y)=k$的个数，即$g(k)={\sum }_{i=1}^n{\sum }_{j=1}^m[gcd(i,j)=k]$
那么重新构造ans的表达式为
$$ans=\sum _{k=1}^{min(n,m)}\sigma (k)g(k)$$

对$g(k)$反演得(关于这个常见反演不会的看这里莫比乌斯反演–懵逼反演系列）

$$g(k)=\sum _{k|d}\mu (\frac{d}{k})\lfloor \frac{n}{d}\rfloor \lfloor \frac{m}{d}\rfloor$$

带入原表达式

$$ans=\sum _{k=1}^{min(n,m)}\sigma (k)\sum _{k|d}\mu (\frac{d}{k})\lfloor \frac{n}{d}\rfloor \lfloor \frac{m}{d}\rfloor$$

更换$d,k$的枚举顺序

$$\sum _{d=1}^{min(n,m)}\lfloor \frac{n}{d}\rfloor \lfloor \frac{m}{d}\rfloor \sum _{k|d}\sigma (k)\mu (\frac{d}{k})$$

现在考虑a的限制，显然只有$\sigma (k)\le a$时其才会产生贡献
所以考虑离线记录询问，按$a$升序排序，预处理所有$\sigma (k)$，升序排序
遍历每个询问时用树状数组动态维护${\sum }_{k|d}\sigma (k)\mu (\frac{d}{k})$这一部分的前缀和

---

#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#include<queue>
using namespace std;
#define lowbit(x) ((x)&(-x))
 
int read()
{
    int f=1,x=0;
    char ss=getchar();
    while(ss<'0'||ss>'9'){if(ss=='-')f=-1;ss=getchar();}
    while(ss>='0'&&ss<='9'){x=x*10+ss-'0';ss=getchar();}
    return f*x;
}
 
const int maxn=100010;
int T;
int miu[maxn];
int prim[maxn],vis[maxn],cnt;
int sum[maxn],ans[maxn];
 
struct node{int id,val; node(){val=0;} }f[maxn];
struct query{int id,n,m,a;}q[maxn];
bool cmp1(query x,query y){ return x.a<y.a;} 
bool cmp2(node x,node y){ return x.val<y.val;}
 
void add(int x,int v){ for(int i=x;i<=maxn-5;i+=lowbit(i)) sum[i]+=v;}
int qsum(int x){ int res=0; for(int i=x;i>0;i-=lowbit(i))res+=sum[i]; return res;}
 
void work(int n)
{
    miu[1]=1;
    for(int i=2;i<=n;++i)
    {
        if(!vis[i]) prim[++cnt]=i,miu[i]=-1;
        for(int j=1;j<=cnt;++j)
        {
            if(i*prim[j]>n) break;
            vis[i*prim[j]]=1;
            if(i%prim[j]==0) break;
            else miu[i*prim[j]]=-miu[i];
        }
    }
     
    for(int i=1;i<=n;++i)
    {
        f[i].id=i;
        for(int j=i;j<=n;j+=i) 
        f[j].val+=i;
    }
    sort(f+1,f+1+n,cmp2);
}
 
int qans(int k)
{
    int n=q[k].n,m=q[k].m;
    int res=0,lim=min(n,m);
    for(int ll=1,rr;ll<=lim;ll=rr+1)
    {
        rr=min(n/(n/ll),m/(m/ll));
        res+=(n/ll)*(m/ll)*(qsum(rr)-qsum(ll-1));
    }
    return res;
}
 
int main() 
{
   T=read(); int lim=0;
   for(int i=1;i<=T;++i)
   {
        q[i].n=read(); q[i].m=read();
        q[i].a=read(); q[i].id=i;
        lim=max(lim,min(q[i].n,q[i].m));
   }
    
   work(lim); int k=1;
   sort(q+1,q+1+T,cmp1); 
   for(int i=1;i<=T;++i)
   {
        while(k<=lim&&f[k].val<=q[i].a)
        {
            for(int j=f[k].id;j<=lim;j+=f[k].id)
            add(j,f[k].val*miu[j/f[k].id]);
            k++;
        }
        ans[q[i].id]=qans(i);
   }
    
   for(int i=1;i<=T;++i)
   printf("%d\n",ans[i]&(~(1<<31)));
   return 0; 
}