[SDOI2014] 数表

Description

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

Solution

先按套路进行一些推导

我们用 BIT 来维护 (h(x))，将所有询问按照 (a) 升序排序，当 (a) 的容许范围扩大时，就将对应的所有 (x) 满足 (f(x)=a) 对 (h) 的贡献进行修改，具体地，对于 (x)，它会对 (ix) 产生 (a mu(i)) 的贡献，这部分单点修改即可。在进行整除分块计算的过程中，在 BIT 上区间求和即可。

#include <bits/stdc++.h>
using namespace std;

#define int long long
const int N = 1000005;
const int MAXN = 1000005;
const int mod = 1ll<<31;

namespace bit {
    const int N = 1000000;
    int ar[N]; // index: 1 ~ N
    int lowbit(int t) { return t & (-t); }
    void add(int i, int v) {
        for (; i < N; ar[i] += v, i += lowbit(i));
    }
    int sum(int i) {
        int s = 0;
        for (; i > 0; s += ar[i], i -= lowbit(i));
        return s;
    }
    int query(int i,int j) {
        return sum(j) - sum(i-1);
    }
}

bool isNotPrime[MAXN + 1];
int mu[MAXN + 1], phi[MAXN + 1], primes[MAXN + 1], cnt, f[N], idx[N];

map<int,int> mp;
vector<int> vec[N];

int n=1e5;

inline void euler() {
    isNotPrime[0] = isNotPrime[1] = true;
    mu[1] = 1;
    for (int i = 2; i <= MAXN; i++) {
        if (!isNotPrime[i]) {
            primes[++cnt] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= cnt; j++) {
            int t = i * primes[j];
            if (t > MAXN) break;
            isNotPrime[t] = true;
            if (i % primes[j] == 0) {
                mu[t] = 0;
                break;
            } else {
                mu[t] = -mu[i];
            }
        }
    }
    for(int i=1;i<=n;i++) {
        int j=1;
        for(;j*j<i;j++) if(i%j==0) {
            f[i]+=j;
            f[i]+=i/j;
        }
        if(j*j==i) {
            f[i]+=j;
        }
        mp[f[i]]++;
    }
    int ind=0;
    for(auto &i:mp) {
        i.second=++ind;
        idx[ind]=i.first;
    }
    for(int i=1;i<=n;i++) {
        vec[mp[f[i]]].push_back(i);
    }
}

struct qry {
    int l,r,a,id;
    bool operator < (const qry & b) {
        return a<b.a;
    }
} qr[N];

void refresh(int x,int a) {
    for(int i=1;i*x<=n;i++) {
        bit::add(i*x,a*mu[i]);
    }
}

void refresh(int a) {
    int t=mp[a];
    for(int i:vec[t]) {
        refresh(i,a);
    }
}

int solve(signed n,signed m) {
    if(n==0 || m==0) return 0;
    int ans=0;
    signed l=1,r=0;
    if(n>m) swap(n,m);
    while(l<=n) {
        r=min(n/(n/l),m/(m/l));
        ans+=bit::query(l,r)*(n/l)*(m/l);
        l=r+1;
    }
    return ans;
}

int res[N];

signed main() {
    euler();
    int t,n,m;
    ios::sync_with_stdio(false);
    cin>>t;
    int cc=0;
    while(t--) {
        ++cc;
        int x,y,z;
        cin>>x>>y>>z;
        qr[cc]={x,y,z,cc};
    }
    sort(qr+1,qr+cc+1);
    int pos=0;
    for(int i=1;i<=cc;i++) {
        auto it=mp.upper_bound(qr[i].a);
        --it;
        int tmp=0;
        if(it!=(--mp.begin())) tmp=it->second;
        while(pos<tmp) {
            ++pos;
            refresh(idx[pos]);
        }
        res[qr[i].id]=solve(qr[i].l,qr[i].r);
    }
    for(int i=1;i<=cc;i++) cout<<(res[i]%mod+mod)%mod<<endl;
}