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

24 4 310 10 5

Sample Output

20148

HINT

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

题解

一道比较恶心的题
我们要求的就是$ans = \sum_{i=1}^{N}\sum_{j=1}^{M}g(gcd(i,j))$，其中g(i)指i的约束和
利用莫比乌斯反演化简得：
$ans = \sum_{T=1}^{N}\lfloor\frac{N}{T}\rfloor\lfloor\frac{M}{T}\rfloor*\sum_{i|T}\mu(\frac{T}{i})g(i)$
然后很常规：
前面部分分块
后面部分维护T的前缀和
维护g(i)的方式：枚举自然数i和i的倍数T，将i的倍数T对应的g(T)加上$\mu(\frac{T}{i})g(i)$
预处理复杂度$O(nlogn)$

但是题目要求我们求<=a<script type="math/tex" id="MathJax-Element-19"><=a</script>的g(i)，我们就将i按照g(i)排序，将询问按照a排序，每次询问前先将前缀和更新到不大于a，此时用树状数组维护前缀和

小技巧：对$2^{31}$取模，可以自然溢出，输出时&上$2^{31}-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;
}