bzoj3529 [Sdoi2014]数表

Description

---

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

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

Solution

---

如果没有a的限制很简单，就是$\sum_{d=1}^{n}{d\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor}$，小范围验证了一下大致是没问题的
考虑有a的限制但是先不管a的限制，那么

$$ans=\sum_{i=1}^{n}\sum_{j=1}^{m}{f(gcd(i,j))}$$ 其中 $f(x)$表示x的因数和
枚举一个gcd，有 $$ans=\sum_{i=1}^{n}{\sum_{j=1}^{m}{\sum_{d=1}^{n}f(d)[gcd(i,j)=d]}}$$
套路一波拉出d来有 $$ans=\sum_{d=1}^{n}{f(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}{\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(i,j)=1]}}$$
一个略微少见的套路强行套一个μ进去 $$ans=\sum_{d=1}^{n}{f(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}{\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{x|gcd(i,j)}{\mu(x)}}}$$
拉出x得到 $$ans=\sum_{d=1}^{n}{f(d)\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\mu(x)\lfloor\frac{n}{dx}\rfloor\lfloor\frac{m}{dx}\rfloor}$$
这里设 $T=dx$，那么又可以枚举T得到 $$ans=\sum_{T=1}^{n}\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum_{d|T}f(d)*\mu(\frac{T}{d})$$
仔细观察发现后面是个狄利克雷卷积，也就是能用线性筛搞出来。考虑加上a的限制，那么就是先求出所有的f并排序，离线询问并按a排序，把所有f[now]<=a的用树状数组插入维护，每次暴力 $O(\sqrt m+\sqrt n)$更新答案即可

推很久错了然后被指出是换∑的时候漏了东西

Code

---

#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <stdlib.h>
#include <time.h>
#define rep(i,st,ed) for (int i=st;i<=ed;++i)
#define lowbit(x) ((x)&(-(x)))

typedef long long LL;
const LL MOD=0x7fffffff;
const int N=200005;
struct Q {int n,m,a,id;} q[N];

LL mu[N+1],prime[N+1],f[N+1];
LL ans[N+1],sum[N+1];
int r[N+1];
bool not_prime[N+1];

void pre_work() {
    mu[1]=1;
    rep(i,2,N) {
        if (!not_prime[i]) {
            prime[++prime[0]]=i;
            mu[i]=-1;
        }
        for (int j=1;j<=prime[0]&&i*prime[j]<=N;j++) {
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0) {
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
    }
    rep(i,1,N) {
        for (int j=i;j<=N;j+=i) {
            f[j]=f[j]+i;
        }
    }
}

void add(int x,LL v) {
    while (x<=N) {
        sum[x]+=v;
        x+=lowbit(x);
    }
}

LL query(int x) {
    LL ret=0;
    while (x) {
        ret+=sum[x];
        x-=lowbit(x);
    }
    return ret;
}

LL solve1(int n,int m) {
    LL ans=0;
    rep(i,1,n) {
        rep(j,1,m) {
            LL tot=0;
            rep(k,1,std:: min(i,j)) {
                if (i%k==0&&j%k==0) tot=(tot+k)%MOD;
            }
            ans=(ans+tot)%MOD;
        }
    }
    return ans;
}

bool cmp1(Q a,Q b) {
    return a.a<b.a;
}

bool cmp2(int x,int y) {
    return f[x]<f[y];
}

int main(void) {
    pre_work();
    int T; scanf("%d",&T);
    rep(i,1,T) {
        scanf("%d%d%d",&q[i].n,&q[i].m,&q[i].a);
        if (q[i].m<q[i].n) std:: swap(q[i].n,q[i].m);
        q[i].id=i;
    }
    rep(i,1,N) r[i]=i;
    std:: sort(q+1,q+T+1,cmp1);
    std:: sort(r+1,r+N+1,cmp2);
    int now=1;
    for (int i=1;i<=T;i++) {
        while (f[r[now]]<=q[i].a) {
            for (int j=r[now];j<=N;j+=r[now]) {
                add(j,f[r[now]]*mu[j/r[now]]);
            }
            now++;
        }
        for (int j=1,k;j<=q[i].n;j=k+1) {
            k=std:: min(q[i].n/(q[i].n/j),q[i].m/(q[i].m/j));
            ans[q[i].id]+=(q[i].n/j)*(q[i].m/j)*(query(k)-query(j-1));
        }
    }
    rep(i,1,T) printf("%lld\n", ans[i]&MOD);
    return 0;
}