BZOJ3529 [Sdoi2014]数表 【莫比乌斯反演】

3529: [Sdoi2014]数表

Time Limit: 10 Sec Memory Limit: 512 MB
Submit: 2151 Solved: 1080
[Submit][Status][Discuss]
Description

有一张N×m的数表，其第i行第j列（1 < =i < =礼，1 < =j < =m）的数值为

能同时整除i和j的所有自然数之和。给定a，计算数表中不大于a的数之和。

Input

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

Output

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

Sample Input

2

4 4 3

10 10 5

Sample Output

20

148

HINT

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

题解

一道比较恶心的题
我们要求的就是ans=∑Ni=1∑Mj=1g(gcd(i,j))，其中g(i)指i的约束和
利用莫比乌斯反演化简得：
ans=∑NT=1⌊NT⌋⌊MT⌋∗∑i|Tμ(Ti)g(i)
然后很常规：
前面部分分块
后面部分维护T的前缀和
维护g(i)的方式：枚举自然数i和i的倍数T，将i的倍数T对应的g(T)加上μ(Ti)g(i)
预处理复杂度O(nlogn)

但是题目要求我们求<=a的g(i)，我们就将i按照g(i)排序，将询问按照a排序，每次询问前先将前缀和更新到不大于a，此时用树状数组维护前缀和

小技巧：对231取模，可以自然溢出，输出时&上231−1【化为正数】

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define LL long long int
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define Redge(u) for (int k = h[u]; k != -1; k = ed[k].nxt)
#define lbt(x) (x & -x)
using namespace std;
const int maxn = 100005,maxm = 20005,INF = 1000000000;
inline int RD(){
    int out = 0,flag = 1; char c = getchar();
    while (c < 48 || c > 57) {if (c == '-') flag = -1; c = getchar();}
    while (c >= 48 && c <= 57) {out = (out << 1) + (out << 3) + c - '0'; c = getchar();}
    return out * flag;
}
int N = 0;
int A[maxn],now = 0,mu[maxn],prime[maxn],primei = 0,Qi,Ans[maxm];
bool isn[maxn];
struct Que{int n,m,a,id;}Q[maxm];
struct Gf{int i,v;}G[maxn];
inline bool operator <(const Que& a,const Que& b){return a.a < b.a;}
inline bool operator <(const Gf& a,const Gf& b){return a.v < b.v;}
inline void add(int u,int v){while (u <= N) A[u] += v,u += lbt(u);}
inline int query(int u){int ans = 0; while (u) ans += A[u],u -= lbt(u); return ans;}
void init(){
    mu[1] = 1;
    for (int i = 2; i <= N; i++){
        if (!isn[i]) prime[++primei] = i,mu[i] = -1;
        for (int j = 1; j <= primei && i * prime[j] <= N; j++){
            isn[i * prime[j]] = true;
            if (i % prime[j] == 0) {mu[i * prime[j]] = 0; break;}
            mu[i * prime[j]] = -mu[i];
        }
    }
    for (int i = 1; i <= N; i++)
        for (int j = i; j <= N; j += i)
            G[j].v += i;
    REP(i,N) G[i].i = i;
    sort(G + 1,G + 1 + N);
}
int main(){
    Qi = RD();
    REP(i,Qi) Q[i].n = RD(),Q[i].m = RD(),Q[i].a = RD(),Q[i].id = i,N = max(N,max(Q[i].n,Q[i].m));
    sort(Q + 1,Q + 1 + Qi);
    init();
    REP(i,Qi){
        while (now < N && G[now + 1].v <= Q[i].a){
            now++;
            for (int j = 1; G[now].i * j <= N; j++)
                add(G[now].i * j,mu[j] * G[now].v);
        }
        int n = Q[i].n,m = Q[i].m; if (n > m) swap(n,m);
        for (int j = 1,nxt; j <= n; j = nxt + 1){
            nxt = min(n / (n / j),m / (m / j));
            Ans[Q[i].id] += (n / j) * (m / j) * (query(nxt) - query(j - 1));
        }
    }
    REP(i,Qi) printf("%d\n",Ans[i] & 0x7fffffff);
    return 0;
}