NKOJ 2841 （SDOI 2014）数表（莫比乌斯反演+树状数组+线性筛）

P2841【SDOI2014 R1D1】数表

问题描述

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

输入格式

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

输出格式

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

样例输入

2
4 4 3
10 10 5

样例输出

20
148

提示

对于30%的数据，1<=n,m<=400,1<=Q<=200
对于另外30%的数据，1<=n,m<=10^5,1<=Q<=10
对于100%的数据，1<=n,m<=10^5,1<=Q<=20000,0<=a<=10^9

---

容易得到对于单次询问

$$Ans=\sum_{i=1}^{n}\sum_{j=1}^{m}\sigma(gcd(i,j)),[\sigma(gcd(i,j))<=a]$$
那么，推导一番 $$Ans=\sum_{d=1}^{N}\sigma(d)\sum_{i=1}^{n}\sum_{j=1}^{m}1[gcd(i,j)=d],N=min(n,m)$$
然后，反演一下 $$Ans=\sum_{d=1}^{N}\sigma(d)\sum_{d|K}^{N}\mu(\frac{K}{d})\lfloor{\frac{n}{K}}\rfloor\lfloor{\frac{m}{K}}\rfloor=\sum_{K=1}^{N}\lfloor{\frac{n}{K}}\rfloor\lfloor{\frac{m}{K}}\rfloor\sum_{d|K}\sigma(d)\mu(\frac{K}{d}),\sigma(d)<=a$$

注意到$\sum_{d|K}\sigma(d)\mu(\frac{K}{d})$是个积性函数，用线性筛预处理，但是只有$\sigma(d)<=a$的项才对答案有贡献，因此离线处理，对询问按照a排序，然后每次将小于等于a的那些项暴力枚举一下倍数，更新到树状数组里面，查询的时候直接用。
复杂度为$O(n\ ln\ n+m\sqrt{n}log_2n)$左右

---

代码：

#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<cstring>
#define ll long long
#define N 100005
using namespace std;
const ll mod=(1ll<<31);
struct node{ll n,m,a,id;}K[N],g[N];
bool cmp(node aa,node bb)
{return aa.a<bb.a;}
ll T,P[N],tot,mu[N],pc[N],pd[N],G[N],Ans[N];
void EU()
{
    ll i,j;mu[1]=1;g[1].a=1;
    for(i=1;i<N;i++)g[i].n=i;
    for(i=2;i<N;i++)
    {
        if(!g[i].a)P[++tot]=i,mu[i]=-1,g[i].a=i+1,pc[i]=i+1,pd[i]=tot;
        for(j=1;j<=tot&&P[j]*i<N;j++)
        if(i%P[j])
        {
            pc[i*P[j]]=P[j]+1;
            pd[i*P[j]]=j;
            mu[i*P[j]]=-mu[i];
            g[i*P[j]].a=g[i].a*pc[i*P[j]]%mod;
        }
        else
        {
            pc[i*P[j]]=pd[i]==j?pc[i]*P[j]+1:P[j]+1;
            pd[i*P[j]]=j;
            g[i*P[j]].a=pd[i]==j?g[i].a/pc[i]*pc[i*P[j]]:g[i].a*pc[i*P[j]];
            pc[i*P[j]]%=mod;
            g[i*P[j]].a%=mod;
        }
    }
    sort(g+1,g+N,cmp);
}
void MD(ll x,ll d)
{for(ll i=x;i<N;i+=(i&-i))G[i]+=d,G[i]%=mod;}
ll GS(ll x)
{
    ll i,sum=0;
    for(i=x;i;i-=(i&-i))sum+=G[i],sum%=mod;
    return sum;
}
int main()
{
    ll i,j,p,k=1,ans;EU();
    scanf("%lld",&T);
    for(i=1;i<=T;i++)scanf("%lld%lld%lld",&K[i].n,&K[i].m,&K[i].a),K[i].id=i;
    sort(K+1,K+T+1,cmp);
    for(i=1;i<=T;i++)
    {
        while(g[k].a<=K[i].a)
        {
            for(j=g[k].n;j<N;j+=g[k].n)MD(j,g[k].a*mu[j/g[k].n]%mod);
            k++;
        }
        ans=0;
        for(j=1;j<=K[i].n&&j<=K[i].m;j=p+1)
        {
            p=min(K[i].n/(K[i].n/j),K[i].m/(K[i].m/j));
            ans+=(K[i].n/j)*(K[i].m/j)%mod*(GS(p)-GS(j-1))%mod;ans%=mod;
        }
        Ans[K[i].id]=(ans+mod)%mod;
    }
    for(i=1;i<=T;i++)printf("%lld\n",Ans[i]);
}