洛谷 P3312 [SDOI2014]数表 莫比乌斯反演+树状数组+线性筛

题目描述

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

输入输出格式

输入格式：
输入包含多组数据。

输入的第一行一个整数Q表示测试点内的数据组数

接下来Q行，每行三个整数n，m，a(|a| < =10^9)描述一组数据。

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

输入输出样例

输入样例#1：
2
4 4 3
10 10 5
输出样例#1：
20
148
说明

1 < =N．m < =10^5 ， 1 < =Q < =2*10^4

分析：i,j的公约数可以看作是gcd(i,j)的约数，设为F(gcd(i,j)),其中F(gcd(i,j))<=a。

$$ans=\sum_{i=1}^n\sum_{j=1}^mF(gcd(i,j))\ mod\ 2^{31}$$
$$=\sum_{D=1}^nF(D)\sum_{d=1}^{\lfloor\frac{n}{D}\rfloor}\mu(d)*\lfloor\frac{n}{Dd}\rfloor*\lfloor\frac{m}{Dd}\rfloor$$

设$T=Dd$,
则

$$\sum_{T=1}^{n}\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum_{D|T}F(D)*\mu(\frac{T}{D})$$

前面显然可以整除分块，关键是后面了。
设

$$g(x)=\sum_{i|x}F(i)*\mu(\frac{x}{i})$$

显然$F=g*id$，就是一个积性函数，线筛一波，顺便把$\mu$也给筛出来。

我们发现$F(i)*\mu(\frac{x}{i})$有贡献时，$F(i)<=a$。

于是我们可以把询问离线，把各个询问的a从小到大排序，也就是说，每次把ai-1到ai的F对g的贡献给算出来加入到g，也就是要个F排序了，那么相当于动态维护小于等于a时的g前缀和，显然这个时候就是树状数组一波了。

代码：

// luogu-judger-enable-o2
#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>

const int maxc=2e4+7;
const int maxn=1e5+7;

using namespace std;

int test,cnt;
int not_prime[maxn],prime[maxn],mul[maxn],low[maxn];
int t[maxn];

struct rec{
    int x,y;
}f[maxn];

struct node{
    int n,m,k,ans,num;
}q[maxc];

bool cmp1(node x,node y)
{
    return x.k<y.k;
}

bool cmp2(rec x,rec y)
{
    return x.y<y.y;
}

bool cmp3(node x,node y)
{
    return x.num<y.num;
}

void getmul(int n)
{
    f[1]=(rec){1,1};
    mul[1]=1;
    for (int i=2;i<=n;i++)
    {
        if (!not_prime[i])
        {
            prime[++cnt]=i;
            f[i]=(rec){i,i+1};
            low[i]=i;
            mul[i]=-1;
        }
        for (int j=1;j<=cnt;j++)
        {
            if (i*prime[j]>n) break;
            not_prime[i*prime[j]]=1;
            if (low[i]%prime[j]==0) low[i*prime[j]]=low[i]*prime[j];
                               else low[i*prime[j]]=prime[j];
            if (i*prime[j]==low[i*prime[j]])
            {
                f[i*prime[j]]=(rec){i*prime[j],f[i].y+i*prime[j]};
                mul[i*prime[j]]=0;
            }
            else
            {
                int x=low[i*prime[j]],y=i*prime[j]/x;
                f[i*prime[j]]=(rec){i*prime[j],f[x].y*f[y].y};
                mul[i*prime[j]]=mul[x]*mul[y];
            }
        }
    }
    sort(f+1,f+n+1,cmp2);
}

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

void ins(int x,int k)
{
    while (x<=1e5)
    {
        t[x]+=k;
        x+=lowbit(x);
    }
}

int getsum(int x)
{
    int ans=0;
    while (x>0)
    {
        ans+=t[x];
        x-=lowbit(x);
    }
    return ans;
}

int main()
{
    scanf("%d",&test);
    getmul(1e5);
    for (int i=1;i<=test;i++)
    {
        scanf("%d%d%d",&q[i].n,&q[i].m,&q[i].k);
        q[i].num=i;
        if (q[i].n>q[i].m) swap(q[i].n,q[i].m);
    }
    sort(q+1,q+test+1,cmp1);
    int d=1;
    for (int i=1;i<=test;i++)
    {
        while ((f[d].y<=q[i].k) && (d<=1e5))
        {
            for (int j=f[d].x;j<=1e5;j+=f[d].x)
            {
                ins(j,f[d].y*mul[j/f[d].x]);
            }
            d++;
        }
        int x=0,y=0;
        int n=q[i].n,m=q[i].m;
        for (int j=1,last;j<=n;j=last+1)
        {
            last=min(n/(n/j),m/(m/j));
            y=getsum(last);
            q[i].ans+=(n/j)*(m/j)*(y-x);
            x=y;
        }
    }
    sort(q+1,q+test+1,cmp3);
    for (int i=1;i<=test;i++) printf("%d\n",q[i].ans&0x7fffffff);
}