Luogu 3312 - [SDOI2014]数表 （莫比乌斯反演、分块、树状数组）

Luogu 3312 - [SDOI2014]数表

题意

给定\(n,m,a\)，求解下式

\[\sum_{i=1}^n\sum_{j=1}^m S(\gcd(i,j))\ [S(\gcd(i,j))\leq a]\ (mod\ 2^{31}) \]

其中\(S(x)\)表示\(x\)的因子和

限制

\(1\leq T\leq 2\times 10^4,\ 1\leq n,m\leq 10^5,\ |a|\leq10^9\)

思路

令\(n\leq m\)

先不考虑\(a\)的限制

将原式转为枚举\(d=\gcd(i,j)\)，得

\[ans=\sum_{x=1}^nS(x)\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=x] \]

引入函数\(g(x)\)

\[g(x)=\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=x] \]

转为枚举\(d=\gcd(i,j)\)，得

\[g(x)=\sum_{x|d}\mu(\frac dx)\lfloor\frac nd\rfloor\lfloor\frac md\rfloor \]

将\(g(x)\)代回原式，得

\[ans=\sum_{x=1}^nS(x)g(x)\\ =\sum_{x=1}^nS(x)\sum_{x|d}\mu(\frac dx)\lfloor\frac nd\rfloor\lfloor\frac md\rfloor \]

明显对于\(S(x)\)是可以预处理出来的，但是后半部分随着\(d\)的变化需要重新计算，再加上询问组数，时间复杂度为\(O(n\times\sqrt n\times m)\)，显然无法得出正解

所以可以转换\(x\)与\(d\)的枚举顺序，将枚举\(x\)的倍数\(d\)转为枚举\(d\)的因子\(x\)，得

\[ans=\sum_{d=1}^n\lfloor\frac nd\rfloor\lfloor\frac md\rfloor\sum_{x|d}S(x)\mu(\frac dx) \]

显然对于\(\lfloor\frac nd\rfloor\lfloor\frac md\rfloor\)可以通过分块来计算

引入函数\(h(x)\)

\[h(x)=\sum_{x|d}S(x)\mu(\frac dx) \]

如果不考虑\(a\)的限制，\(h(x)\)可以直接通过预处理前缀和获得，那么本题就算完成了

但考虑到\(a\)的限制后，我们可以先记录所有的询问，按照\(a\)从小到大的顺序去处理

由于\(a\)有序，我们可以在每次询问处理完毕后，转到下一个询问时，把\(a\)增加的部分所满足的\(h(x)\)加入前缀和中，再进行分块

为了使得修改前缀和更为高效，所以这里需要使用树状数组来维护

如果处理到询问\(i\)时，就将所有\(a_{i-1}\lt S(x)\leq a_i\)的\(S(x)\)找出，将\(S(x)\mu(i)\)加入树状数组的\(ix,\ 2ix,\ 3ix...\)这些位置内，再配合分块即可

为了使查找的功能更为高效，预处理完\(S\)数组后带编号排序一遍即可

最后，模数比较特殊，直接使用带符号\(32\)位整型（\(int\)）自然溢出即可

如果最后算出的答案自然溢出为负数，加上\(2^{31}\)即可（注意使用\(long\ long\)类型加）

时间复杂度：预处理\(O(n\sqrt n)\)，树状数组维护\(O(nlogn)\)，分块处理多组数据\(O(m\sqrt n)\)，总时间复杂度为\(O((n+m)\sqrt n+nlogn)\)

代码

Case Max (912ms/1500ms)

O2 Case Max (391ms/1500ms)

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=100000,M=20000;

int primcnt;
ll mu[N+50],prim[N+50];
bool vis[N+50];
ll ans[M+50];

struct S_func
{
    int S,id;
    bool operator < (const S_func &hs) const
    {
        return S<hs.S;
    }
}S[N+50];

struct que
{
    int n,m,a,id;
    bool operator < (const que &nd) const
    {
        return a<nd.a;
    }
}q[M+50];

struct BIT
{
    typedef int Type;
    Type a[N+50];
    void upd(int p,Type v){for(;p<=N;p+=(p&-p))a[p]+=v;}
    Type qry(int p){Type r=0;for(;p;p-=(p&-p))r+=a[p];return r;}
    Type qry(int l,int r){return qry(r)-qry(l-1);}
}tree;

void init(int n)
{
    memset(tree.a,0,sizeof tree.a);
    memset(vis,false,sizeof vis);
    mu[1]=1;
    primcnt=0;
    for(int j=1;j<=n;j++)
    {
        S[j].id=j;
        S[j].S=1;
    }
    for(int i=2;i<=n;i++)
    {
        if(!vis[i])
        {
            prim[primcnt++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<primcnt;j++)
        {
            int k=i*prim[j];
            if(k>n)
                break;
            vis[k]=true;
            if(i%prim[j]==0)
            {
                mu[k]=0;
                break;
            }
            mu[k]=-mu[i];
        }
        
        for(int j=i;j<=n;j+=i)
            S[j].S+=i;
    }
    sort(S+1,S+1+n);
}

int main()
{
    init(N);
    
    int T;
    scanf("%d",&T);
    for(int t=1;t<=T;t++)
    {
        scanf("%d%d%d",&q[t].n,&q[t].m,&q[t].a);
        q[t].id=t;
    }
    sort(q+1,q+1+T);
    
    for(int t=1,p=1;t<=T;t++)
    {
        int &n=q[t].n,&m=q[t].m,&a=q[t].a;
        if(n>m)
            swap(n,m);
        
        while(p<=N&&S[p].S<=a) //满足条件的都加入BIT
        {
            for(int i=S[p].id;i<=N;i+=S[p].id)
                tree.upd(i,S[p].S*mu[i/S[p].id]);
            p++;
        }
        
        int res=0;
        for(int x=1;x<=n;)
        {
            int nxt=min(n/(n/x),m/(m/x));
            res=res+(n/x)*(m/x)*tree.qry(x,nxt);
            x=nxt+1;
        }
        if(res<0)
            res+=(1LL<<31); //自然溢出结果若为负数记得加个模数
        ans[q[t].id]=res;
    }
    
    for(int t=1;t<=T;t++)
        printf("%d\n",ans[t]);
    
    return 0;
}