题目地址：洛谷：【P3312】[SDOI2014]数表 – 洛谷、BZOJ：Problem 3529. — [Sdoi2014]数表

题目描述

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

输入输出格式

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

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

输入输出样例

输入样例#1：

2
4 4 3
10 10 5

输出样例#1：

20
148

说明

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

题解

参考资料：题解 P3312 【[SDOI2014]数表】 – ldlsn 的博客 – 洛谷博客
神奇的反演。
冷静分析一下，如果固定i和j，如何求能同时整除i和j的所有自然数之和。规定F(x)为x的因数和，例如F(2)=1+2=3, F(6)=1+2+3+6=12，则上面这个问题的答案就是F[\mathrm{gcd}(i, j)]。
得到了上面的结论后，不考虑a的限制，我们可以知道答案就是\sum_{i=1}^n \sum_{j=1}^m F[\mathrm{gcd}(i, j)]，下面就是对这个答案式的变形。
我们知道常用的反演结论g(x) = \sum_{i=1}^n \sum_{j=1}^m [\mathrm{gcd}(i, j) == x] = \sum_{i|d} \mu \left( \frac{d}{i} \right) \left\lfloor \frac{n}{d} \right\rfloor \left\lfloor \frac{m}{d} \right\rfloor，直接套进答案里
\sum_{i=1}^n F(i)g(i) = \sum_{i=1}^n F(i) \sum_{i|d} \mu \left( \frac{d}{i} \right) \left\lfloor \frac{n}{d} \right\rfloor \left\lfloor \frac{m}{d} \right\rfloor
改变一下求和顺序
\sum_{d=1}^n \left\lfloor \frac{n}{d} \right\rfloor \left\lfloor \frac{m}{d} \right\rfloor \sum_{i|d} F(i) \mu \left( \frac{d}{i} \right)
如果我们可以计算出f(x) = \sum_{i|x} F(i) \mu(\frac{x}{i})的值前缀和，那么可以除法分块做，而F(x)可以枚举倍数求和，这些都不是问题。
接下来我们来看带上限制应该怎么做，考虑离线处理，把询问按a限制从小到大排序，每一次将符合a条件的F(x)对f(x)产生的贡献扔进一个维护前缀和的数据结构中（比如树状数组），为回答询问的除法分块求前缀和做准备。
至于那个奇怪的模数，自然溢出后&0x7fffffff就好。
总复杂度是O(n \sqrt{n} \log n)。

代码

// Code by KSkun, 2018/6
#include <cstdio>
#include <cctype>

#include <algorithm>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) 
        ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1;
    register char c = fgc();
    while(!isdigit(c)) {
        if(c == '-') neg = -1;
        c = fgc();
    }
    while(isdigit(c)) {
        res = (res << 1) + (res << 3) + c - '0';
        c = fgc();
    }
    return res * neg;
}

const int MAXN = 100005;

int prime[MAXN], mu[MAXN], ptot;
bool notprime[MAXN];

inline void sieve() {
    notprime[1] = true; mu[1] = 1;
    for(int i = 2; i < MAXN; i++) {
        if(!notprime[i]) {
            prime[++ptot] = i;
            mu[i] = -1;
        }
        for(int j = 1; j <= ptot && i * prime[j] < MAXN; j++) {
            int t = i * prime[j];
            notprime[t] = true;
            if(i % prime[j]) {
                mu[t] = -mu[i];
            } else {
                mu[t] = 0; break;
            }
        }
    }
}

int bit[MAXN];

inline int lowbit(int x) {
    return x & -x;
}

inline void add(int x, int v) {
    for(int i = x; i < MAXN; i += lowbit(i)) {
        bit[i] += v;
    }
}

inline int query(int x) {
    int res = 0;
    for(int i = x; i; i -= lowbit(i)) {
        res += bit[i];
    }
    return res;
}

struct Node {
    int x, v;
    inline bool operator<(const Node &rhs) const {
        return v < rhs.v;
    }
} F[MAXN];

inline void calf() {
    for(int i = 1; i <= 100000; i++) {
        F[i].x = i;
        for(int j = 1; i * j <= 100000; j++) {
            F[i * j].v += i;
        }
    }
    std::sort(F + 1, F + 100001);
}

inline int calc(int n, int m) {
    int ans = 0;
    if(n > m) std::swap(n, m);
    for(int i = 1; i <= n;) {
        int j = std::min(n / (n / i), m / (m / i));
        ans += (query(j) - query(i - 1)) * (n / i) * (m / i);
        i = j + 1;
    }
    return ans;
}

struct Query {
    int n, m, a, id;
} qur[MAXN];

inline bool qcmp(Query a, Query b) {
    return a.a < b.a;
}

int T, ans[MAXN];

int main() {
    sieve(); calf();
    T = readint();
    for(int i = 1; i <= T; i++) {
        qur[i].n = readint(); qur[i].m = readint(); qur[i].a = readint(); qur[i].id = i;
    }
    std::sort(qur + 1, qur + T + 1, qcmp);
    int now = 0;
    for(int i = 1; i <= T; i++) {
        while(now < 100000 && F[now + 1].v <= qur[i].a) {
            now++;
            for(int j = 1; j * F[now].x <= 100000; j++) {
                add(j * F[now].x, mu[j] * F[now].v);
            }
        }
        ans[qur[i].id] = calc(qur[i].n, qur[i].m);
    }
    for(int i = 1; i <= T; i++) {
        printf("%d\n", ans[i] & 0x7fffffff);
    }
    return 0;
}