bzoj3529: [Sdoi2014]数表 莫比乌斯反演

bzoj3529: [Sdoi2014]数表

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

分析

首先，我们令

$$f(i,j)=\sum_{d|i\bigwedge d|j}d$$
则 $$ans=\sum_i^n\sum_i^m[f(i,j)\leq a]f(i,j)$$
考虑f，不难发现 $$f(i.j)=\sum_{d|gcd(i,j)}d=\sigma(gcd(i,j))$$
$其中\sigma(i)为i的约束和$
首先咱们先不考虑a的限制
$$ans=\sum_i^n\sum_i^m\sigma(gcd(i,j))$$
然后枚举约数
$$ans=\sum_d^n\sigma(d)\sum_i^n\sum_i^m [gcd(i,j))==d]$$
根据模型(详见： bzoj2820YY的GCD)得
$$\sum_i^n\sum_i^m [gcd(i,j))==d]=\sum_k^{\lfloor\frac{n}{d}\rfloor}\mu(k)\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor$$
于是有
$$ans=\sum_d^n\sigma(d)\sum_k^{\lfloor\frac{n}{d}\rfloor}\mu(k)\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor$$
转换枚举变量，令 $D=kd$
$$ans=\sum_D^n\lfloor\frac{n}{D}\rfloor\lfloor\frac{m}{D}\rfloor \sum_{d|D}\sigma(d)\mu(\frac{D}{d})$$
这个时候，我们发现我们还没有考虑a的限制。。
我们考虑最开始，发现a本身的限制就是加在f(i,j)也就是 $\sigma(d)$上的，所以我们最后的答案可以直接写成
$$ans=\sum_D^n\lfloor\frac{n}{D}\rfloor\lfloor\frac{m}{D}\rfloor \sum_{d|D}\sigma(d)\mu(\frac{D}{d})[\sigma(d) \leq a]$$
那么，对于 $\sum \limits_D^n\lfloor\frac{n} {D}\rfloor\lfloor\frac{m}{D}\rfloor$可以使用下底函数分块。
而对于 $\sum \limits_{d|D}\sigma(d)\mu(\frac{D}{d})[\sigma(d) \leq a]$ 处理条件 $[\sigma(d) \leq a]$，我们离线将a排序，从小到大把 $\sigma(d)\mu(\frac{D}{d})$的插入树状数组即可。使用树状数组的原因

一是随着a的改变这玩意的值也会改变，需要支持单点修改。二是分块求答案的时候需要做一个区间求和。……时间复杂度$O(T\sqrt n log_n)$。——form xxcc

代码

具体细节看代码，本题还是有很多值得推敲之处的。

#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const int N = 110000;
const int M = 440000;
const int inf = 0x7fffffff;
int read() {
    char ch = getchar(); int x = 0, f = 1;
    while(ch < '0' || ch > '9') {if(ch == '-') f = -1; ch = getchar();}
    while(ch >= '0' && ch <= '9') {x = x * 10 + ch - '0'; ch = getchar();}
    return x * f;
}
pair<int, int>f[N];
int mu[N], p[N], t[N], a[N], id[N], qn[N], qm[N], ans[N], mx, tot, top, m;
bool vis[N];
void add(int x, int add) {for(int i = x; i <= mx; i += i & -i) t[i] += add;}
int query(int x) {int ret = 0; for(int i = x; i; i -= i & -i) ret += t[i]; return ret;}
bool cmp(int b, int c) {return a[b] < a[c];}

void mobius() {
    mu[1] = 1;
    for(int i = 2;i <= mx; ++i) {
        if(!vis[i]) {p[++tot] = i; mu[i] = -1;}
        for(int j = 1;j <= tot && i * p[j] <= mx; ++j) {
            vis[i * p[j]] = true;
            if(i % p[j]) mu[i * p[j]] = -mu[i];
            else break;
        }
    } 
    for(int i = 1;i <= mx; ++i)
        for(int j = i;j <= mx; j += i)
            f[j].first += i;
    for(int i = 1;i <= mx; ++i) f[i].second = i;
}

void work(int id) {
    while(top <= mx && f[top].first <= a[id]) {
        for(int i = f[top].second; i <= mx; i += f[top].second)
            add(i, f[top].first * mu[i / f[top].second]);
        ++top;
    }
    int n = qn[id], m = qm[id];
    for(int i = 1, pos;i <= n; i = pos + 1) {
        pos = min(n / (n / i), m / (m / i));
        ans[id] += (n / i) * (m / i) * (query(pos) - query(i - 1));
    }

}

int main() {
    m = read();
    for(int i = 1;i <= m; ++i) {
        qn[i] = read(); qm[i] = read(); if(qn[i] > qm[i]) swap(qn[i], qm[i]);
        a[i] = read(); id[i] = i; mx = max(mx, qn[i]);
    }
    mobius();
    sort(id + 1, id + m + 1, cmp);
    sort(f + 1, f + mx + 1); top = 1;
    for(int i = 1;i <= m; ++i) work(id[i]);
    for(int i = 1;i <= m; ++i) printf("%d\n", ans[i] & inf);
    return 0;
}