📚【模板】数论分块

简介

用于以 \(O(\sqrt n)\) 求解类似于 \(\sum\limits_{i = 1}^{\min_{j = 1}^{k}\{n_j\}}f(i)\times\prod\limits_{j = 1}^{k}\left\lfloor\frac{n_j}{i}\right\rfloor\) 的式子。

前提是可以方便的求出 \(f(i)\) 的前缀和。

#define llt long long int
llt NTS(llt* pre,llt upd,llt p) {
	llt res = 0LL;
	for(register llt l = 1, r;l <= upd;++l) {
		r = upd/(upd/l);
		res = (res+(upd/l)*(pre[r]-pre[l-1]+p)%p)%p;
	}
	return res;
}

例题

\(\textrm{luogu P2260 [清华集训2012]模积和}\)

求解 \(\sum\limits_{i = 1}^{n}\sum\limits_{j = 1}^{m}(n \bmod i)\cdot(m \bmod j)\;,\;(i \not= j)\)

本来不大难，但是 \(i\not=j\) 有点小麻烦。

\[\begin{aligned}ans &= \sum\limits_{i = 1}^{n}\left(n-\left\lfloor\frac{n}{i}\right\rfloor\cdot i\right)\times\sum\limits_{j = 1}^{m}\left(m-\left\lfloor\frac{m}{j}\right\rfloor\cdot j\right)-\sum\limits_{i = 1}^{n}\left(n-\left\lfloor\frac{n}{i}\right\rfloor\cdot i\right)\cdot\left(m-\left\lfloor\frac{m}{i}\right\rfloor\cdot i\right)\\ &= \left(n^2-\sum\limits_{i = 1}^{n}\left\lfloor\frac{n}{i}\right\rfloor\cdot i\right)\times\left(m^2-\sum\limits_{j = 1}^{m}\left\lfloor\frac{m}{j}\right\rfloor\cdot j\right)-\sum\limits_{i = 1}^{n}\left(n\cdot m-n\cdot \left\lfloor\frac{m}{i}\right\rfloor\cdot i-m\cdot \left\lfloor\frac{n}{i}\right\rfloor\cdot i+\left\lfloor\frac{n}{i}\right\rfloor\cdot \left\lfloor\frac{m}{i}\right\rfloor\right)\end{aligned}\]

数论分块即可。

顺便 \(\sum\limits_{i = 1}^{n}i = \frac{n\times(n+1)}{2}\) ， \(\sum\limits_{i = 1}^{n}i^2 = \frac{n\times(n+1)\times(2\cdot n+1)}{6}\) 。

#include <iostream>
#define llt long long int
const llt moyn = 19940417;
const llt inv_2 = 9970209;
const llt inv_6 = 3323403;
llt n, m;
llt l, r;
llt res_1, res_2, res_3;
llt tmp_1, tmp_2, tmp_3;
llt pre_1(llt upd) {
	return upd*(upd+1LL)%moyn*inv_2%moyn;
}
llt pre_2(llt upd) {
	return upd*(upd+1LL)%moyn*((upd<<1)+1LL)%moyn*inv_6%moyn;
}
llt bmod(llt x) {
	return (x%moyn+moyn)%moyn;
}
int main() {
	scanf("%lld %lld",&n,&m);
	res_1 = n*n%moyn;
	res_2 = m*m%moyn;
	for(l = 1;l <= n;l = r+1) {
		r = n/(n/l);
		res_1 = bmod(res_1-(n/l)*bmod(pre_1(r)-pre_1(l-1)));
	}
	for(l = 1;l <= m;l = r+1) {
		r = m/(m/l);
		res_2 = bmod(res_2-(m/l)*bmod(pre_1(r)-pre_1(l-1)));
	}
	for(l = 1;l <= std :: min(n,m);l = r+1) {
		r = std :: min(n/(n/l),m/(m/l));
		tmp_1 = (r-l+1LL)*n%moyn*m%moyn;
		tmp_2 = ((n/l)*m+(m/l)*n)%moyn*bmod(pre_1(r)-pre_1(l-1))%moyn;
		tmp_3 = (n/l)*(m/l)%moyn*bmod(pre_2(r)-pre_2(l-1))%moyn;
		res_3 = bmod(res_3+tmp_1-tmp_2+tmp_3);
	}
	printf("%lld\n",bmod(res_1*res_2-res_3));
	return 0;
}