[bzoj3529][Sdoi2014]数表

Description

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

Input

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

Output

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

Sample Input

2
4 4 3
10 10 5

Sample Output

20
148

Hint

1 < =N．m < =10^5​5​​ ， 1 < =Q < =2×10410^410​4

在水了很多很多的莫比乌斯反演水题之后还是这个玩意让我停了下来

嗯我做到的第一个离线的数论题，想了一阵呢

因为莫比乌斯反演水的题挺多的了，套路都一样，具体的化简方法和yy的gcd一样，就不再细讲了

根据题意，表格中的每个数a[i][j]都表示gcd(i,j)的约数和

预处理f(i)表示i的约数和，复杂度O（n log n），也可以线筛时弄出来，但是这题的数据范围不大，前者更好实现些

那么我们要求的式子就是：Σ_(i=1->n) Σ_(j=1->m) f(gcd(i,j)) [f(gcd(i,j))<=a]

改变枚举顺序：Σ_{(p=1->min(n,m))} f(p) Σ_(i=1->n/p) Σ_(j=1->m/p) [gcd(i,j)==p] [f(p)<=a]

反演：Σ_{(p=1->min(n,m))} f(p) Σ_{(d=1->min(n,m)/p)} μ(d)*(n/d/p)*(m/d/p)

化简改变枚举项：Σ_{(dp=1->min(m,n))} (n/dp)*(m/dp) Σ_(p|dp) f(p)*μ(dp/p)

显然对于同一组询问，画下划线的部分的每一个dp值可以预处理，设之为g(dp)

那么结合数论分块，处理g(dp)的前缀和，就可以O(4√n)回答每一组询问

接下来考虑每一组询问中不同的a值如何处理。

当f函数的值超过了a的时候，显然可以直接将其f值变为0，不会对答案产生影响。

我们尝试将询问按a值从大到小排序，然后预处理出所有的f和g值。

 接下来开始处理每个询问前，找到所有超过a的f值，删除其对g值的影响，即g(i*j)-=f(i)*μ(j)

因为a值已经被排序，每个数的f值都会被至多一次删除影响，次数为(n log n)

而为了提高找到f值的效率，我们可以在预处理时将f排序。

下一个问题是我们需要查询g值的区间和，需要的操作是：区间查询，单点修改。

可以使用树状数组来存g值，因为其自带前缀和性质，可以满足我们的需要。

共(n log n)次修改，(T*4√n*2)次查询，每次的复杂度都是O(log n)

总复杂度O(n log²n+8T log n √n)，可以接受。

另一种类似但更优的做法是询问按从小到大排序，那样就不用预处理所有g了。

在每次询问前不是删除影响而是产生影响，这么做能节约常数时间并且无需预处理。

A之后才想出来，我好笨啊。

附第一种想法的笨重代码：

#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
const long long mod=2147483648ll;
long long ans[20005],f[100005],g[100005],tr[100005];short miu[100005],np[100005];int p[100005],nop,sqr,t,m,n,a,qt=100000,res;
void add(int p,long long num){for(;p<=100005;p+=p&-p)tr[p]=(tr[p]+num+mod)%mod;}
long long sum(int p,long long ans=0){for(;p;p-=p&-p)ans=(ans+tr[p]+mod)%mod;return ans;}
inline int min(int a,int b){return a<b?a:b;}
struct qs{
    int m,n,a,num;
    friend bool operator<(qs z,qs b){
        return z.a>b.a;
    }
}q[20005];
struct fff{
    int ff,num;
    friend bool operator<(fff a,fff b){
        return a.ff<b.ff;
    }
}fs[100005];
int main(){
    miu[1]=1;
    for(int i=2;i<=100000;++i){
        if(!np[i])p[++nop]=i,miu[i]=-1;
        for(int j=1;j<=nop&&p[j]*i<=100000;++j)
            if(i%p[j])np[i*p[j]]=1,miu[i*p[j]]=miu[i]*-1;
            else{np[i*p[j]]=1;break;}
    }
    for(int i=1;i<=100000;++i){sqr=sqrt(i);
        for(int j=1;j<=sqr;++j)if(i%j==0)f[i]=f[i]+j+i/j;
        if(sqr*sqr==i)f[i]=f[i]-sqr;
    }
    for(int i=1;i<=100000;++i)fs[i].ff=f[i],fs[i].num=i;
    sort(fs+1,fs+100001);qt=100000;
    for(int i=1;i<=100000;++i)for(int j=1;i*j<=100000;++j)g[i*j]=(g[i*j]+f[i]*miu[j]+mod)%mod;
    for(int j=1;j<=100000;++j)add(j,g[j]);
    scanf("%d",&t);
    for(int i=1;i<=t;++i)scanf("%d%d%d",&q[i].m,&q[i].n,&q[i].a),q[i].num=i;
    sort(q+1,q+1+t);q[0].a=100000;
    for(int i=1;i<=t;++i){
        if(q[i].a!=q[i-1].a)while(fs[qt].ff>q[i].a){
            for(int j=1;j*fs[qt].num<=100000;++j)add(j*fs[qt].num,-fs[qt].ff*miu[j]);
            qt--;
        }
        res=q[i].m<q[i].n?q[i].m:q[i].n;
        for(int ii=1,last;ii<=res;ii=last+1)last=min(q[i].m/(q[i].m/ii),q[i].n/(q[i].n/ii)),ans[q[i].num]=(ans[q[i].num]+(q[i].m/ii)*(q[i].n/ii)*(sum(last)-sum(ii-1))%mod)%mod;
        ans[q[i].num]=(ans[q[i].num]+mod)%mod;
    }
    for(int i=1;i<=t;++i)printf("%lld
",ans[i]);
}