Luogu 3312 [SDOI2014]数表

在这一篇里把所有的套路写全方便自己之后复习。

首先是一个小学生数学：$a$整除$b$ $ = $  $frac{b}{a}$

也就是说这题中格子$(i, j)$的值就是既能被$i$整除又能被$j$整除的所有自然数的和。

小学数学没学好……

那么思考一下就能列出本题的式子：

　　　　　　$sum_{i = 1}^{n}sum_{j = 1}^{m}f(gcd(i, j))$       假设$n leq m$并且$f(gcd(i, j)) leq a$

$f(i)$表示$i$的约数和

发现$a$的限制十分讨厌，先不考虑它把式子推出来。

先枚举$d$

　　　　　　$sum_{d = 1}^{n}f(d)sum_{i = 1}^{n}sum_{j = 1}^{m}[gcd(i, j) == d]$

后面直接枚举$d$的倍数

　　　　　　$sum_{d = 1}^{n}f(d)sum_{i = 1}^{left lfloor frac{n}{d} ight floor}sum_{j = 1}^{left lfloor frac{m}{d} ight floor}[gcd(i, j) == 1]$

在这里可以把$[gcd(i, j) == 1]$进行单位元代换了

　　　　　　$sum_{d = 1}^{n}f(d)sum_{i = 1}^{left lfloor frac{n}{d} ight floor}sum_{j = 1}^{left lfloor frac{m}{d} ight floor}sum_{t | gcd(i, j)}mu (t)$

枚举$t$

　　　　　　$sum_{d = 1}^{n}f(d)sum_{t = 1}^{left lfloor frac{n}{d} ight floor}mu (t)sum_{i = 1}^{left lfloor frac{n}{d} ight floor}sum_{j = 1}^{left lfloor frac{m}{d} ight floor}[t|gcd(i, j)]$

后面直接枚举$t$的倍数

　　　　　　$sum_{d = 1}^{n}f(d)sum_{t = 1}^{left lfloor frac{n}{d} ight floor}mu (t)left lfloor frac{n}{td} ight floorleft lfloor frac{m}{td} ight floor$

枚举$T = dt$

　　　　　　$sum_{T = 1}^{n}left lfloor frac{n}{T} ight floorleft lfloor frac{m}{T} ight floorsum_{d | T}mu (frac{T}{d})f(d)$

记$h(T) = sum_{d | T}mu (frac{T}{d})f(d)$

　　　　　　$sum_{T = 1}^{n}left lfloor frac{n}{T} ight floorleft lfloor frac{m}{T} ight floor h(T)$

$h(i)$显然可以线性筛，如果不考虑$a$的情况下已经可以整除分块$sqrt{n}$求解了。

现在考虑$a$的限制，发现只有当$f(i) leq a$的时候才会对答案产生贡献，我们可以将所有询问离线，按照$a$从小到大排序，然后写一个支持单点修改区间查询的数据结构来做这个产生贡献的$h(i)$的前缀和的计算。

当然是常数小又好写的$bit$啦。

时间复杂度$O(qsqrt{n}log MaxN + MaxNlogMaxN)$。

在本题中取模只要让$int$自然溢出，最后再& $2^{31} - 1$。

Code：

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef pair <int, int> pin;

const int N = 1e5 + 5;
const int M = 2e4 + 5;

int qn, maxn = 0, pCnt = 0, pri[N], low[N], mu[N], d[N];
bool np[N];
pin f[N];

struct Querys {
    int n, m, r, id, ans;
} q[M];

bool cmp(const Querys &x, const Querys &y) {
    return x.r < y.r;
}

inline void read(int &X) {
    X = 0; char ch = 0; int op = 1;
    for(; ch > '9' || ch < '0'; ch = getchar())
        if(ch == '-') op = -1;
    for(; ch >= '0' && ch <= '9'; ch = getchar())
        X = (X << 3) + (X << 1) + ch - 48;
    X *= op;
}

inline void chkMax(int &x, int y) {
    if(y > x) x = y;
}

inline int min(int x, int y) {
    return x > y ? y : x;
} 

inline void sieve() {
    mu[1] = d[1] = 1;
    for(int i = 2; i <= maxn; i++) {
        if(!np[i]) pri[++pCnt] = i, mu[i] = -1, d[i] = 1 + i, low[i] = i; 
        for(int j = 1; j <= pCnt && pri[j] * i <= maxn; j++) {
            np[i * pri[j]] = 1;
            if(i % pri[j] == 0) {
                low[i * pri[j]] = low[i] * pri[j];
                mu[i * pri[j]] = 0;
                if(low[i] == i) d[i * pri[j]] = d[i] * pri[j] + 1;
                else d[i * pri[j]] = d[i / low[i]] * d[pri[j] * low[i]];
                break;
            }
            low[i * pri[j]] = pri[j];
            mu[i * pri[j]] = -mu[i];
            d[i * pri[j]] = d[i] * d[pri[j]];
        }
    }
    
/*    for(int i = 1; i <= 20; i++)
        printf("%d ", mu[i]);
    printf("
");
    for(int i = 1; i <= 20; i++)
        printf("%d ", d[i]);
    printf("
");    */
    
    for(int i = 1; i <= maxn; i++) 
        f[i].first = d[i], f[i].second = i;
    sort(f + 1, f + 1 + maxn);
}

namespace BinaryIndexTree {
    int s[N];
    
    #define lowbit(p) (p & (-p))
    
    inline void modify(int p, int v) {
        for(; p <= maxn; p += lowbit(p))
            s[p] += v;
    }
    
    inline int query(int p) {
        int res = 0;
        for(; p > 0; p -= lowbit(p))
            res += s[p];
        return res;
    } 
    
} using namespace BinaryIndexTree;

inline void solve() {
    sort(q + 1, q + 1 + qn, cmp);
    int now = 0; 
    for(int i = 1; i <= qn; i++) {
        int rep = min(q[i].n, q[i].m), res = 0;
        
        for(; now + 1 <= maxn && f[now + 1].first <= q[i].r; ) {
            ++now;
            for(int j = f[now].second; j <= maxn; j += f[now].second)
                modify(j, f[now].first * mu[j / f[now].second]);
        }
        
        for(int l = 1, r; l <= rep; l = r + 1) {
            r = min((q[i].n / (q[i].n / l)), (q[i].m / (q[i].m / l)));
            res += (q[i].n / l) * (q[i].m / l) * (query(r) - query(l - 1));
        }
        q[q[i].id].ans = res;
    }
}

int main() {
    read(qn);
    for(int i = 1; i <= qn; i++) {
        read(q[i].n), read(q[i].m), read(q[i].r);
        q[i].id = i;
        chkMax(maxn, q[i].n), chkMax(maxn, q[i].m);
    }        
    
    sieve();    
    solve();
    
    for(int i = 1; i <= qn; i++)
        printf("%d
", q[i].ans & (0x7fffffff));
    return 0;
}