[SDOI2014]数表

给反演题写题解好累啊。。。

因为给反演题写题解都要打很多latex，所以虽然之前做了一些反演题，但都懒得写题解。借这篇题解写一下反演的套路。

传送门

设$f(i)$为i的约数和，求$\sum_{i=1}^n\sum_{j=1}^m[f(gcd(i,j)<=a)]f(gcd(i,j))$（默认n<=m）

不考虑a
$$
\sum_{i=1}^n\sum_{j=1}^mf(gcd(i,j))=\sum_{i=1}^n\sum_{j=1}^m\sum_{d=1}^n[gcd(i,j)=d]f(d)
$$
把d放到前面，再枚举d的倍数
$$
\sum_{d=1}^nf(d)\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)==1]
$$
枚举i,j的因子k，反演，再把k提到前面，改成枚举k的倍数
$$
\sum_{d=1}^{n}f(d)\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}\sum_{k|(i,j)}μ(k)
=\sum_{d=1}^{n}f(d)\sum_{k=1}^{n/d}μ(k)\sum_{i=1}^{n/dk}\sum_{j=1}^{m/dk}
$$
枚举$d\times k=T$
$$
\sum_{T=1}^{n}\sum_{d|T}f(d)\timesμ(\frac{T}{d})\times \frac{n}{T}\times \frac{m}{T}
$$
之后$\frac{n}{T} \times \frac{m}{T}$可以整除分块，$g(T)=\sum_{d|T}f(d)\timesμ(\frac{T}{d})$是个关于T的函数，线筛出f和μ，埃筛求g即可。

线性筛约数和的方法

那么到这里你已经解决了一个普通的反演题，我也把这篇题解最麻烦的部分写完了。

对于a的问题，其实是让大于a的f(x)不参与g(x)的更新即可，所以做法如下：

把询问离线，按a排序，每次把目前小于a的f(x)加到x的所有倍数的g函数中，用树状数组维护g的前缀和。

每个f(x)最多只会被加入一次，最坏复杂度是$O(n\ log^2n+q\sqrt n\ log\ n)$。

#include<cstdio>
#include<iostream>
#include<algorithm>
using namespace std;
struct node{
    int n,m,a,h;
}q[20001];
struct delta{
    int z,h;
}f[200001];
int o[100001],sum[1000001],low[100001];
int prime[100001],t,ans[20001],maxn,mu[100001];
inline bool cmp(node c,node d){return c.a<d.a;}
inline bool cmp2(delta c,delta d){return c.z<d.z;}
inline int lowbit(int x){return x&(-x);}
inline void add(int x,int y){while (x<=maxn) sum[x]+=y,x+=lowbit(x);}
inline int query(int x){int num=0; while (x) num+=sum[x],x-=lowbit(x); return num;}
inline void euler(){
    o[1]=f[1].z=mu[1]=low[1]=1;
    for (int i=2; i<=maxn; i++){
        if (!o[i]){
            prime[++prime[0]]=i;
            f[i].z=low[i]=i+1; ; mu[i]=-1;
        }
        for (int j=1; j<=prime[0]&&prime[j]*i<=maxn; j++){
            o[i*prime[j]]=1;
            if (i%prime[j]==0){
                f[i*prime[j]].z=f[i].z/low[i]*(low[i]*prime[j]+1);
                low[i*prime[j]]=low[i]*prime[j]+1;
                break;
            }
            f[i*prime[j]].z=f[i].z*f[prime[j]].z;
            low[i*prime[j]]=prime[j]+1;
            mu[i*prime[j]]=-mu[i];
        }
    }
    for (int i=1; i<=maxn; i++) f[i].h=i;
//	for (int i=1; i<=maxn; i++){
//		g[i].h=i;
//		for (int j=i; j<=maxn; j+=i)
//			g[j].z+=f[i].z*mu[j/i];
//	}
//	for (int i=1; i<=maxn; i++) printf("%d ",g[i].z);
//	puts("");
//	sort(g+1,g+maxn+1,cmp2);
    sort(f+1,f+maxn+1,cmp2);
}
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);
        if (q[i].n>q[i].m) swap(q[i].n,q[i].m);
        q[i].h=i; maxn=max(maxn,q[i].m);
    }
    sort(q+1,q+t+1,cmp); euler(); int k=1,n,m,a;
    for (int i=1; i<=t; i++){
        n=q[i].n,m=q[i].m,a=q[i].a;
        while (f[k].z<=a&&k<=maxn){
            for (int j=f[k].h; j<=maxn; j+=f[k].h)
                add(j,f[k].z*mu[j/f[k].h]);
//				g[j]+=;
//			add(g[k].h,g[k].z);
            k++;
        }
        for (int l=1,r; l<=n; l=r+1){
            r=min(n/(n/l),m/(m/l));
            ans[q[i].h]+=(n/l)*(m/l)*(query(r)-query(l-1));
        }
    }
    for (int i=1; i<=t; i++) printf("%d\n",ans[i]&(~(1<<31)));
    return 0;
}