[莫反] [数据结构] P3312 [SDOI2014] 数表

posted on 2024-03-14 08:54:03 | under | source

先不考虑 \(\le a\) 的限制，莫反可得：\(ans=\sum\limits_{T=1}^{\min(n,m)}\frac nT\frac mTH(T)\)。其中 \(H(T)=\sum\limits_{d|T} \sigma(d)\mu(\frac Td)\)。

现在加上限制，将 \(\le a\) 融入式子看起来不太可行，这启示我们从函数值入手。如果能快速算出 \(\le a\) 时 \(\sigma(d)\) 的值（\(d>a\) 就令其值为 \(0\)），然后直接整除分块就好了。

暴力肯定是不行的。如果离线做，使得 \(a,d\) 有序，那么显然一个 \(d\) 至多被加入一次。

现在考虑加入一个 \(d\) 的影响。会对所有 \(d\mid T\) 的 \(H(T)\) 造成 \(\sigma(d)\mu(\frac Td)\) 的影响，直接修改即可。由于整除分块时算的是前缀和，所以把 \(H\) 放在树状数组上维护即可。

众所周知 \(\sum\limits_{i=1}^n \frac ni\approx n\log n\)，所以维护 \(H\) 复杂度 \(O(q\log n^2)\)，分块总复杂度 \(O(q\sqrt n \log n)\)。

代码

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

#define int long long
const int N = 1e5 + 5, mod = 1ll << 31;
int q, o[N], mu[N], prime[N], pc, hs[N], H[N], ans[N];
bool flag[N];
int idd[N], idw[N], t[N];
struct ask{int n, m, a, id;} b[N];
struct sigma{int res, id;} s[N];

inline bool cmp(ask A, ask B) {return A.a < B.a;}
inline bool cmp2(sigma A, sigma B) {return A.res < B.res;}
inline void upd(int a, int k) {k = (k % mod + mod) % mod; for(; a < N; a += a & -a) t[a] = (t[a] + k) % mod;} //树状数组维护H前缀和 
inline int qry(int a) {int res = 0; for(; a > 0; a -= a & -a) res = (res + t[a]) % mod; return res;}
inline void init(){
	o[1] = mu[1] = flag[1] = 1;
	for(int i = 2; i < N; ++i){
		if(!flag[i]) {o[i] = 1 + i, mu[i] = -1, prime[++pc] = i, hs[i] = i;}
		for(int j = 1; j <= pc && i * prime[j] < N; ++j){
			int x = i * prime[j];
			mu[x] = i % prime[j] ? -mu[i] : 0;
			hs[x] = i % prime[j] ? prime[j] : hs[i] * prime[j];
			o[x] = i % prime[j] ? o[i] * (1 + prime[j]) : o[i] + hs[x] * o[x / hs[x]];
			flag[x] = 1;
			if(i % prime[j] == 0) break;
		}
	}
	for(int i = 1; i < N; ++i) s[i] = {o[i], i};
	sort(s + 1, s + N, cmp2);
}
signed main(){
	init();
	cin >> q;
	for(int i = 1; i <= q; ++i) scanf("%lld%lld%lld", &b[i].n, &b[i].m, &b[i].a), b[i].id = i;
	sort(b + 1, b + 1 + q, cmp);
	for(int i = 1, now = 1; i <= q; ++i){
		while(now < N && s[now].res <= b[i].a){
			for(int j = s[now].id; j < N; j += s[now].id) upd(j, mu[j / s[now].id] * s[now].res);
			++now; 
		}
		if(b[i].n > b[i].m)	swap(b[i].n, b[i].m);
		for(int l = 1; l <= b[i].n;){
			int r = min(b[i].n / (b[i].n / l), b[i].m / (b[i].m / l));
			ans[b[i].id] = (ans[b[i].id] + (b[i].n / l) * (b[i].m / l) % mod * ((qry(r) - qry(l - 1)) + mod) % mod) % mod;
			l = r + 1;
		}
	}
	for(int i = 1; i <= q; ++i) printf("%lld\n", ans[i]);
	return 0;
}