LibreOJ 2193 「SDOI2014」数表

显然题目所求为 n∑i=1m∑j=1[σ(gcd(i,j))≤a]σ(gcd(i,j))。
先不考虑 a 的限制，推一推式子。

n∑i=1m∑j=1σ(gcd(i,j))=min(n,m)∑i=1σ(i)n∑x=1m∑y=1[gcd(x,y)=i]=min(n,m)∑i=1σ(i)∑d|iμ(id)⌊nd⌋⌊md⌋=min(n,m)∑i=1⌊ni⌋⌊mi⌋∑d|iσ(d)μ(id)

现在我们来考虑 a 的限制。
设 sum(i)=∑d|iσ(d)μ(id)，那么显然对于任意的 σ(i)≤a，它只会对满足 i|d 的 sum(d) 有贡献。
于是可以把 σ(i) 排序，并把询问离线按 a 排序，然后随着 a 的增大，不可能会有 σ(i) 的贡献消失，只可能会增加，所以枚举倍数并更新贡献即可。
这样的话，增加贡献的次数是 O(n∑i=1ni)≈O(nlogn)。

接着考虑如何维护贡献。
由于要单点修改并且区间查询，必须在线，常数要求较高，当然首选 BIT。

此外，因为模数很特殊，所以自然溢出并在最后按位与 231−1 即可。

好奇为什么约数和这种东西都要线性筛（
反正 O(nlogn) 随便跑，懒得推了。

代码：

#include <cstdio>
#include <algorithm>
#define lowbit(x) ((x) & -(x))
using namespace std;
const int N = 1e5;
const int Q = 2e4;
int q;
struct note
{
    int n,m,a,id;
    inline bool operator<(const note &o) const
    {
        return a < o.a || (a == o.a && id < o.id);
    }
} a[N + Q + 5];
int cnt,prime[N + 5],vis[N + 5],mu[N + 5];
long long c[N + 5],ans[Q + 5];
inline void update(int x,long long k)
{
    for(;x <= N;x += lowbit(x))
        c[x] += k;
}
inline long long query(int x)
{
    long long ret = 0;
    for(;x;x -= lowbit(x))
        ret += c[x];
    return ret;
}
int main()
{
    mu[1] = 1,a[1].id = -1;
    for(register int i = 1;i <= N;++i)
    {
        a[i].id = -i;
        for(register int j = 1;i * j <= N;++j)
            a[i * j].a += i;
    }
    for(register int i = 2;i <= N;++i)
    {
        if(!vis[i])
            prime[++cnt] = i,mu[i] = -1;
        for(register int j = 1;j <= cnt && i * prime[j] <= N;++j)
        {
            vis[i * prime[j]] = 1;
            if(!(i % prime[j]))
                break;
            mu[i * prime[j]] = -mu[i];
        }
    }
    scanf("%d",&q);
    for(register int i = 1;i <= q;++i)
        scanf("%d%d%d",&a[N + i].n,&a[N + i].m,&a[N + i].a),a[N + i].id = i;
    sort(a + 1,a + N + q + 1);
    for(register int i = 1,tot = 0;i <= N + q;++i)
        if(a[i].id < 0)
        {
            a[i].id = -a[i].id;
            for(register int j = 1;a[i].id * j <= N;++j)
                update(a[i].id * j,(long long)a[i].a * mu[j]);
        }
        else
        {
            if(a[i].n > a[i].m)
                swap(a[i].n,a[i].m);
            for(register int l = 1,r;l <= a[i].n;l = r + 1)
            {
                r = min(a[i].n / (a[i].n / l),a[i].m / (a[i].m / l));
                ans[a[i].id] += (query(r) - query(l - 1)) * (a[i].n / l) * (a[i].m / l);
            }
            if(++tot == q)
                break;
        }
    for(register int i = 1;i <= q;++i)
        printf("%lld\n",ans[i] & (1LL << 31) - 1);
}