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

测试地址：数表
做法：本题需要用到莫比乌斯反演+树状数组。
首先忽略$a$的限制，我们要求的是：
$ans=\sum_{i=1}^n\sum_{j=1}^mD(\gcd(i,j))$
其中$D(i)$指$i$的约数和。
我们把上式改成枚举公因数$d$，不妨设$n\le m$，则：
$ans=\sum_{d=1}^nD(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]$，我们发现这就是求$\gcd$等于$d$的数对个数（这简直是句废话）。
又令$F(d)=\sum_{d|i}f(i)$，显然$F(d)$为$\gcd$等于$d$倍数的数对个数，那么也有$F(d)=\lfloor \frac{n}{d}\rfloor\lfloor \frac{m}{d}\rfloor$。根据莫比乌斯反演定理的第二种形式，有：
$f(d)=\sum_{d|i}\mu(\frac{i}{d})F(i)=\sum_{d|i}\mu(\frac{i}{d})\lfloor \frac{n}{i}\rfloor\lfloor \frac{m}{i}\rfloor$
把这个式子带进$ans$那个式子，有：
$ans=\sum_{d=1}^nD(d)\sum_{d|i}\mu(\frac{i}{d})\lfloor \frac{n}{i}\rfloor\lfloor \frac{m}{i}\rfloor$
交换$d,i$的位置，有：
$ans=\sum_{i=1}^n\lfloor \frac{n}{i}\rfloor\lfloor \frac{m}{i}\rfloor\sum_{d|i}D(d)\mu(\frac{i}{d})$
显然我们如果预处理出$g(i)=\sum_{d|i}D(d)\mu(\frac{i}{d})$的前缀和，就可以用数论分块处理每个询问了。
那么现在我们考虑$a$的限制，实际上$a$是在限制只有$D(i)\le a$的$D(i)$对$g(i)$有贡献，因此我们把所有询问按$a$从小到大排序，对新产生的贡献暴力在$g$中单点修改，因为$D(i)$对$g(j)$产生贡献当且仅当$i|j$，那么修改次数显然是$n\log n$级别的，然后要在数论分块中维护$g$的前缀和查询，这个显然可以用常数小又好写的树状数组解决。
那么我们就解决了此题，时间复杂度为$O(Q\sqrt n\log n+n\log^2n)$。注意到这题模数很特殊，直接用int自然溢出，最后输出时对$2^{31}-1$取个按位与就行了。
以下是本人代码：

#include <bits/stdc++.h>
using namespace std;
int T,p[100010];
int maxn,D[100010],mu[100010],s[100010],ans[20010];
int prime[100010];
bool vis[100010]={0};
struct query
{
    int id,n,m,a;
}q[20010];

bool cmp(query a,query b)
{
    return a.a<b.a;
}

bool cmpp(int a,int b)
{
    return D[a]<D[b];
}

int lowbit(int x)
{
    return x&(-x);
}

void add(int x,int c)
{
    for(int i=x;i<=maxn;i+=lowbit(i))
        s[i]+=c;
}

int sum(int x)
{
    int ans=0;
    for(int i=x;i;i-=lowbit(i))
        ans+=s[i];
    return ans;
}

void init()
{
    scanf("%d",&T);
    maxn=0;
    for(int i=1;i<=T;i++)
    {
        scanf("%lld%lld%lld",&q[i].n,&q[i].m,&q[i].a);
        if (q[i].n>q[i].m) swap(q[i].n,q[i].m);
        maxn=max(maxn,q[i].n);
        q[i].id=i;
    }
    sort(q+1,q+T+1,cmp);

    for(int i=1;i<=maxn;i++)
        for(int j=1;i*j<=maxn;j++)
            D[i*j]=D[i*j]+i;
    for(int i=1;i<=maxn;i++)
        p[i]=i;
    sort(p+1,p+maxn+1,cmpp);
}

void calc_mu()
{
    mu[1]=1;
    prime[0]=0;
    for(int i=2;i<=maxn;i++)
    {
        if (!vis[i])
        {
            prime[++prime[0]]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=prime[0]&&i*prime[j]<=maxn;j++)
        {
            vis[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
    }
}

int main()
{
    init();
    calc_mu();

    int nowp=1;
    for(int i=1;i<=T;i++)
    {
        while(D[p[nowp]]<=q[i].a)
        {
            for(int j=1;j*p[nowp]<=maxn;j++)
                add(j*p[nowp],D[p[nowp]]*mu[j]);
            nowp++;
        }
        int n=q[i].n,m=q[i].m;
        ans[q[i].id]=0;
        for(int j=n;j>=1;j=max(n/(n/j+1),m/(m/j+1)))
        {
            int l=max(n/(n/j+1)+1,m/(m/j+1)+1),r=j;
            ans[q[i].id]=ans[q[i].id]+(sum(r)-sum(l-1))*(n/j)*(m/j);
        }
    }

    for(int i=1;i<=T;i++)
        printf("%d\n",ans[i]&2147483647);

    return 0;
}