題目

有一張 n×m 的數表，其第 i 行第 j 列（1 <= i <= n, 1 <= j <= m）的數值為
能同時整除 i 和 j 的所有自然數之和。給定 a , 計算數表中不大於 a 的數之和。

題解

先不考慮不大於a這一限制條件。

令\(f(x)\)為x的約數和，那麼數表(i,j)的值即為\(f(\gcd(i, j))\)。

開始的處理思路類似YY的gcd

\[\sum\limits^n_{i=1}{\sum\limits^m_{j=1}{f(\gcd(i, j))}} \]

\[\sum\limits^{\min(n,m)}_{d_1}{\sum\limits^{\frac{n}{d_1}}_{i=1}{\sum\limits^{\frac{m}{d_1}}_{j=1}{[gcd(i,j)=1]f(d_1)}}} \]

\[\sum\limits^{\min(n,m)}_{d_1}{\sum\limits^{\frac{n}{d_1}}_{i=1}{\sum\limits^{\frac{m}{d_1}}_{j=1}{\sum\limits_{d_2|\gcd(i,j)}{\mu(d_2)}f(d_1)}}} \]

\[\sum\limits^{\min(n,m)}_{d_1}{\sum\limits_{d_2}{\mu(d_2)f(d_1) \lfloor\frac{n}{pd}\rfloor\lfloor\frac{m}{pd}\rfloor}} \]

令\(T=d_1d_2\)

\[\sum\limits^{min(n, m)}_{T=1}{\sum\limits_{d|T}{\mu(\frac{T}{d})f(d)\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor}} \]

\[\sum\limits^{min(n, m)}_{T=1}{\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum\limits_{d|T}{\mu(\frac{T}{d})f(d)}} \]

令\(F(T)=\sum\limits_{d|T}{\mu(\frac{T}{d})f(d)}\)。注意，式子中的f(d)若大於a則為零。如果我們能快速求出在限制下，\(F(T)\)的字首和，那麼就可以以\(O(\sqrt{n})\)的複雜度處理每個詢問。

這裡將詢問按a從小到大排序，將f(d)從小到大排序。若當前f(i)小於等於a，則將其貢獻加入（將其加入每個i的倍數的F中）。由於a的值單調，一共加n次，每次加\(\frac{n}{i}\)

，複雜度\(O(nln(n))\)。然後用樹狀陣列維護字首和。時間複雜度\(O(nln(n)+Q*\sqrt{n})\)。

#include <iostream>
#include <map>
#include <algorithm>

#define endl '\n'
#define IOS std::ios::sync_with_stdio(0); cin.tie(0); cout.tie(0)
#define FILE freopen(".//data_generator//in.txt","r",stdin),freopen("res.txt","w",stdout)
#define FI freopen(".//data_generator//in.txt","r",stdin)
#define FO freopen("res.txt","w",stdout)
#define mp make_pair
#define seteps(N) fixed << setprecision(N) 
typedef long long ll;

using namespace std;
/*-----------------------------------------------------------------*/

ll gcd(ll a, ll b) {return b ? gcd(b, a % b) : a;}
#define INF 0x3f3f3f3f

const int N = 3e5 + 10;
const double eps = 1e-5;
const ll M = (1ll<<31) - 1;
bool vis[N];
ll prime[N];
ll cnt;
ll g[N];   
ll sum[N];
struct node {
    ll v;
    int p;
    bool operator<(const node&rhs) const {
        return v < rhs.v;
    }
} f[N];
node ask[N];
node nm[N];
ll tarr[N];
int mu[N];
ll ans[N];



void init() {
    mu[1] = sum[1] = 1;
    for(ll i = 2; i < N; i++) {
        if(!vis[i]) {
            mu[i] = -1;
            g[i] = i + 1;
            sum[i] = g[i];
            vis[i] = 1;
            prime[cnt++] = i;
        }
        for(ll j = 0; j < cnt && prime[j] * i < N; j++) {
            ll p = prime[j];
            vis[p * i] = 1;
            if(i % p == 0) {
                g[p * i] = g[i] * p + 1;
                sum[p * i] = sum[i] / g[i] * g[p * i];
                mu[p * i] = 0;
                break;
            } else {
                sum[p * i] = sum[p] * sum[i];
                g[p * i] = p + 1;
                mu[p * i] = -mu[i];
            }
        }
    }
    
    for(int i = 1; i < N; i++) {
        f[i].p = i;
        f[i].v = sum[i];
    }
    sort(f + 1, f + N);
}

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

void add(int p, ll val) {
    while(p < N) {
        tarr[p] += val;
        p += lowbit(p);
    }
}

ll get(int p) {
    ll res = 0;
    while(p) {
        res += tarr[p];
        p -= lowbit(p);
    }
    return res;
}

int main() {
    IOS;
    init();
    int q;
    cin >> q;
    for(int i = 1; i <= q; i++) {
        ll n, m, a;
        cin >> n >> m >> a;
        ask[i] = node{a, i};
        nm[i] = node{n, m};
    }
    sort(ask + 1, ask + 1 + q);
    int p = 1;
    for(int i = 1; i <= q; i++) {
        while(p < N && f[p].v <= ask[i].v) {
            for(int i = f[p].p; i < N; i += f[p].p) {
                add(i, f[p].v * mu[i / f[p].p]);
            }
            p++;
        }
        ll res = 0;
        int cur = 1;
        int n = nm[ask[i].p].v, m = nm[ask[i].p].p;
        while(cur <= min(n, m)) {
            int p1 = n / (n / cur);
            int p2 = m / (m / cur);
            int nt = min(min(p1, p2), min(n, m));
            res += 1ll * (get(nt) - get(cur - 1)) * (n / cur) * (m / cur);
            cur = nt + 1;
        }
        ans[ask[i].p] = (res & M);
    }
    for(int i = 1; i <= q; i++) cout << ans[i] << endl;
}