【BZOJ3529】【SDOI2014】数表 （莫比乌斯反演+树状数组）

传送门

Description

有一张$n imes m$的数表，其第$i$行第$j$列 $(1≤i≤n,1≤j≤m)$ 的数值为能同时整除$i$和$j$的所有自然数之和。
现在给定$a$，计算数表中不大于$a$的数之和。

Input

输入包含多组数据。
输入的第一行一个整数$Q$表示测试点内的数据组数，接下来Q行，每行三个整数$n,m,a(a≤109)$描述一组数据。

Output

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

题解：

我数学太水了！！又是一道推公式的题：

egin{aligned}
ans&=sum^n_{i=1}sum^m_{j=1}sum_{dmid gdc(i,j)}d\
&=sum_d^nsum^{lfloorfrac{n}{d} floor}_{i=1}sum^{lfloorfrac{m}{d} floor}_{j=1}f(d)[gcd(i,j)=1] &(f(d)=sum_{kmid d}k)\
&=sum_d^n f(d)sum^{lfloorfrac{n}{d} floor}_{i=1}sum^{lfloorfrac{m}{d} floor}_{j=1}sum_{kmid gcd(i,j)}mu(k)\
&=sum_k^nsum_d^n f(d)mu(k)lfloorfrac{n}{kd} floorlfloorfrac{m}{kd} floor\
&设 T=kd \
ans&=sum_T^nsum_{dmid T}f(d)mu(frac{T}{d})lfloorfrac{n}{T} floorlfloorfrac{m}{T} floor\
end{aligned}

前面同样是数论分块，后面依旧是线性筛。

考虑f怎么求

有个约数和定理

若$n=sum_{i=1}^k p_i^{a_i}$，$p_i$为$n$的质因数，那么$n$的约数和 $f(n)$满足$f(n)=prod_{i=1}^ksum_{j=0}^{a_i}p_i^j$

$f$的话首先是个积性函数，我们在筛$mu$的时候想顺便把这个也筛出来

考虑$f(d)$的值，如果说d是质数的话答案显然是$d+1$，下面讨论$d$为合数的情况

设$d=i*p$，其中$p$为质数

$p mid i$，那么$p$和$i$互质，所以$f(d)=f(p)*f(i)$

$p∣i$，设$i=t*p^x$ ，那么根据约数和定理，我们可以得
$$f(i*p) = f(t)f(p^{x+1}) = f(t)sumlimits_{i=0}^{x+1}p^i$$

然后我们把p⁰（也就是1）拿出来，得到
$$f(i * p) = f(t) + f(t)*sumlimits_{i=1}^{x+1}p^i = f(t) + f(t)*f(p^x)*p$$

然后$i = t * p^x$，所以$f(t) * f(p^x) = f(i)$

所以最后就是$f(d) = f(t) + f(i) * p$

然后就可以筛出$f(d)$啦

剩下的东西

现在加上a的限制，其实就是离线处理

我们先将所有的询问按照a的大小排序，然后从小到大处理

因为分块的时候我们要用到的是g(x)的前缀和，所以用一个树状数组来处理

将f(x)排个序，枚举的时候只枚举到f(x)<a，然后枚举另一个约数求出g，丢到树状数组里面去

求答案的时候直接查询就好了

 CODE：

#include<iostream>
#include<cstdio>
#include<algorithm>
using namespace std;

#define lowbit(x) (x&-x)
int T,cnt=0,maxn,ans[20005],tmp[100005];
int f[100005],pri[100005],c[100005],mu[100005];
bool vis[100005];
struct Que{
    int n,m,a,id;
    bool operator<(const Que &b)const{
        return a<b.a;
    }
}q[20005];

bool comp(int a,int b){return f[a]<f[b];}

void init(){
    mu[1]=f[1]=1;
    for(int i=2;i<=maxn;i++){
        if(!vis[i]){
            f[i]=i+1,mu[i]=-1;
            pri[++cnt]=i;
        }
        for(int j=1;j<=cnt&&i*pri[j]<=maxn;j++){
            vis[i*pri[j]]=true;
            if(i%pri[j]==0){
                mu[i*pri[j]]=0;
                int t=i;
                while(t%pri[j]==0)t/=pri[j];
                f[i*pri[j]]=f[t]+pri[j]*f[i];
                break;
            }else{
                f[i*pri[j]]=f[i]*f[pri[j]];
                mu[i*pri[j]]=-mu[i];
            }
        }
    }
}

void add(int x,int y){
    while(x<=maxn){
        c[x]+=y;
        x+=lowbit(x);
    }
}

int sum(int x){
    int ans=0;
    while(x>=1){
        ans+=c[x];
        x-=lowbit(x);
    }
    return ans;
}

int main(){
    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;
        if(q[i].n>q[i].m)swap(q[i].n,q[i].m);
        maxn=max(maxn,q[i].n);
    }
    init();
    sort(q+1,q+T+1);
    for(int i=1;i<=maxn;i++)tmp[i]=i;
    sort(tmp+1,tmp+maxn+1,comp);
    for(int i=1,now=1;i<=T;i++){
        int n=q[i].n,m=q[i].m,a=q[i].a;
        while(now<=maxn&&f[tmp[now]]<=a){
            for(int k=1;k*tmp[now]<=maxn;k++)
                add(k*tmp[now],f[tmp[now]]*mu[k]);
            now++;
        }
        for(int j=1,pos;j<=n;j=pos+1){
            pos=min(n/(n/j),m/(m/j));
            ans[q[i].id]+=(n/j)*(m/j)*(sum(pos)-sum(j-1));
        }
    }
    for(int i=1;i<=T;i++)
        printf("%d
",ans[i]&0x7fffffff);
}