数学提高模板

数学提高

此部分由我前队友完成。浏览体验可能会有所不同因为水平比我强


有些地方无法正常浏览，可能是因为我队友有些用的LateX语法，markdown渲染不出来

常用数学技巧和公式

1. 常用公式

完全立方公式

$(a+b)^3=a^3+3a^2b+3ab^2+b^3=a^3+b^3+3ab(a+b)$

$a^3+b^3=(a+b)^3-3ab(a+b)=(a+b)(a^2+b^2-ab)$

$(a-b)^3=a^3-3a^2b+3ab^2-b^3=a^3-b^3+3ab(b-a)$

$a^3-b^3=(a-b)^3+3ab(a-b)=(a-b)(a^2+b^2+ab)$

（当然，使用二项式定理也能自然的推出来，但是在使用数学公式构造时可能没那么快想起来对应的公式）

1. 拆数

$888..888=8\times 111..111=8\times \frac{999..999}{9}=8\times \frac{10^x-1}{9}$

（也就是相同的某个数构成的数字，可以分解成$111...\times k$的形式，再把$111...$变成$\frac{10^x-1} 9$的形式。

$2334..899=111..111+1111+1111+...$

（也就是一个从左往右的每一位不降的数，可以被分解成若干个$111..$相加

$(x-1)^3+(x-1)^3-x^3-x^3=6x$

102780J的神秘构造公式，具体见后面列出的题目。

任何一个大于2的奇数$x$都可以写成完全平方差应该叫这个吧大概的形式：$x=a^2-b^2$

假设$x=2n+1$，显然$2n+1=(n+1)^2-n^2$

随机化

mt19937生成随机数

在评测机上，$RAND\_MAX$很小，只有32767，$rand()$和$rand\_shuffle()$的函数尽量不要使用。而$Mersenne~Twister$算法提供了另一种随机数生成方法，用到的质数为$2^{19937}-1$，这也是名字的由来。在随机化题目中，真的会有毒瘤出题人卡掉普通的$rand()$随机数算法。

函数用法：

// 获得一个unsigned int类型的随机数
mt19937 rng(time(0));
printf("%u\n", rng());
// 代替random_shuffle()
shuffle(a + 1, a + n + 1, rnd);

完整代码

#include <cstdio>
#include <chrono>
#include <random>
using namespace std;

mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());

int main(){
    printf("%u\n", rng());
    return 0;
}

矩阵乘法

矩阵的封装，重载运算符

例题：p1306 斐波那契数列 计算$\gcd(f_n, f_m)$。

根据数学基础篇的斐波那契数列那一节的性质，$\gcd(f_n,f_m)=f(\gcd(n,m))$，直接计算，这里就可以顺手练一下矩阵封装的模板了（虽然没什么必要）

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define sz 2
const ll N = 2, M = 1e8;

ll read(){
    ll x = 0, f=1;char ch;
    do{ch = getchar();if (ch == '-') f=-1;}while(ch<'0' || ch>'9');
    do{x = x*10 + (ch-'0');ch = getchar();}while(ch>='0' && ch<='9');
    return x*f;
}

struct mat{
	ll a[sz][sz];
	
	inline mat() { memset(a, 0, sizeof(a)); }
	inline mat(ll _a[sz][sz]) {
		for (int i=0;i<sz;++i) {
			for (int j=0;j<sz;++j) {
				a[i][j] = _a[i][j];
			}
		}
	}
	inline void operator = (mat x) {
		for (int i=0; i<sz;++i) {
			for (int j=0;j<sz;++j) {
				a[i][j] = x.a[i][j];
			}
		}
	}
	inline mat operator + (const mat& T) const {
		mat res;
		for (int i=0;i<sz;++i){
			for (int j=0;j<sz;++j){
				res.a[i][j] = (a[i][j] + T.a[i][j]) % M;
			}
		}
		return res;
	}
	inline mat operator - (const mat& T) const {
		mat res;
		for (int i=0;i<sz;++i){
			for (int j=0;j<sz;++j){
				res.a[i][j] = (a[i][j] - T.a[i][j]) % M;
			}
		}
		return res;
	}
	// 只能处理两个矩阵都为sz相同的方阵的乘法!!
	inline mat operator * (const mat& T) const {
		mat res; ll r;
		for (int i=0;i<sz;++i){
			for (int k=0;k<sz;++k){
				r = a[i][k];
				for (int j=0;j<sz;++j){
					res.a[i][j] = (res.a[i][j] + T.a[k][j] * r % M) % M;
				}
			}
		}
		return res;
	}
	inline mat operator ^ (ll x) const {
		mat res, b;
		for (int i=0;i<sz;++i) res.a[i][i] = 1;
		for (int i=0;i<sz;++i){
			for (int j=0;j<sz;++j){
				b.a[i][j] = a[i][j] % M;
			}
		}
		while (x){
			if (x & 1) res = res * b;
			b = b * b;
			x >>= 1;
		}
		return res;
	}
};

// 计算c = a * b
void mul(ll c[], ll a[], ll b[][N]){
	ll t[N] = {0};
	for (int i=0;i<N;++i){
		for (int j=0;j<N;++j){
			t[i] = (a[j] * b[j][i] % M + t[i]) % M;
		}
	}
	memcpy(c, t, sizeof(t));
}

ll Fibonacci(ll n){
	ll _a[N][N] = {
		{0, 1},
		{1, 1}
	};
	mat A(_a);
	ll f1[N] = {1, 1};
	A = A ^ (n - 1);
	mul(f1, f1, A.a);
	return f1[0];
}

ll gcd(ll a, ll b){
	return b == 0 ? a : gcd(b, a % b);
}

int main(void){
	ll n, m;
	n = read(), m = read();
	ll GCD = gcd(n, m);
	printf("%lld\n", Fibonacci(GCD));
	
	return 0; 
}

不封装 直接写函数计算$\times \div a^k$

void mul(int c[], int a[], int b[][N]){
	int t[N] = {0};
	for (int i=0;i<N;++i){
		for (int j=0;j<N;++j){
			t[i] = (1LL * a[j] * b[j][i] + t[i]) % m;
		}
	}
	memcpy(c, t, sizeof(t));
}

void mul(int c[][N], int a[][N], int b[][N]){
	int t[N][N] = {0};
	for (int i=0;i<N;++i){
		for (int j=0;j<N;++j){
			for (int k=0;k<N;++k){
				t[i][j] = (1LL * a[i][k] * b[k][j] + t[i][j]) % m; 
			}
		}
	}
	memcpy(c, t, sizeof(t));
}

void qmi(int a[][N], int b, int p){
	int t[N][N];
	for (int i=0;i<N;++i){
		for (int j=0;j<N;++j){
			t[i][j] = (i == j);
		}
	}
	while (b){
		if (b & 1) mul(t, t, a);
		mul(a, a, a);
		b >>= 1;
	}
	memcpy(a, t, sizeof(t));
}

矩阵乘法应用

(1) 前置知识：矩阵快速幂

和快速幂一样，只不过对象是矩阵。

(2) 应用范围：递推式操作

一般表现为如下形式：某个行向量$X_n=[a_n~~ a_{n-1}~..a_{n-k+1}]$，边界$X_1=[a_{k-1}~~a_{k-2}~..a_0]$，可以找到某个矩阵$A$，满足：$X_n=X_{n-1}\times A$，则可以得到$X_n=X_1\times A^{n-1}$。注意：一定要保证$A$中没有变量。

由于矩阵具有结合律，我们可以先求出$B=A^{n-1}\mod P$，然后再计算$X_1\times B$，即可求出$X_n$，$X_n$的第一个元素就是$a_n$。

(3) 时间复杂度：$O(\log n)$ 但是大常数ORZ。

(4) 例题+模板

① AcWing1303 斐波那契前n项和

现有Fibonacci数列如下：$f_1=1,f_2=1,f_n=f_{n-2}+f_{n-1}$。

定义$S_n=\sum\limits _{i=1}^n f_i$，对于给定的$n,m$，求$S_n \mod m$。

数据范围：$1\leq n \leq 2\times 10^9;1\leq m \leq 10^9 + 10$

首先啊，我们知道对于任意$n>2$，$f_n$只和它的前两项有关，并且，$S_{n+1}=S_{n}+f_{n+1}$，

那么只要我们知道了$f_{n},f_{n+1},S_n$，就可以推出$f_{n+2},S_{n+1}$的值了。

因此定义$X_n=[f_n,f_{n+1},S_n]$，则$X_{n+1}=[f_{n+1},f_{n+2},S_{n+1}]$

那么$[f_n,f_{n+1},S_n]\times \begin{bmatrix} 0~~1~~0 \\ 1~~1~~1 \\ 0~~0~~1 \end{bmatrix}=[f_{n+1},f_{n+2},S_{n+1}]$

因此，$A= \begin{bmatrix} 0~~1~~0 \\ 1~~1~~1 \\ 0~~0~~1 \end{bmatrix}$，求出$A^{n-1}\mod m$后，计算$X_1\times A^{n-1}\mod m$即可。

int n, m;
const int N = 3;

int main(void){
	scanf("%d%d", &n, &m);
	int f1[N] = {1, 1, 1};
	int a[N][N] = {
		{0, 1, 0}, 
		{1, 1, 1}, 
		{0, 0, 1}
	};
	qmi(a, n - 1, m);
	mul(f1, f1, a);
	
	printf("%d\n", f1[2]);
	
	return 0;
}

② 佳佳的斐波那契

这次是求解$T(n)=(f_1+2f_2+..+nf_n)\mod m$的值，数据范围：$1\leq n,m\leq 2^{31}-1$

当所求$T(n)$不能被直接构造时，可以考虑用两个可以算出的数值构造出$T(n)$的答案。

由上题可以看出，$f_n,f_{n+1},S_n$都是可以比较轻松算出来的。

又$\begin{cases} nS_n-T_n=(n-1)f_1+(n-2)f_2+..+f_{n-1} \\ (n+1)S_{n+1}-T_{n+1}=nf_1+(n-1)f_2+..+2f_{n-1}+f_n \end{cases}$

定义$P_n=nS_n-T_n$可得$P_{n+1}-P_n=S_n$，而$T_n=nS_n-P_n$。

下面就很好算了，定义$X_n=[f_n,f_{n+1},S_n,P_n]$，$X_{n+1}=[f_{n+1},f_{n+2},S_{n+1},P_{n+1}]$

则$[f_n,f_{n+1},S_n,P_n]\times \begin{bmatrix} 0~~1~~0~~0\\ 1~~1~~1~~0\\ 0~~0~~1~~1\\ 0~~0~~0~~1 \end{bmatrix}=[f_{n+1},f_{n+2},S_{n+1},P_{n+1}]$

因此，$A=\begin{bmatrix} 0~~1~~0~~0\\ 1~~1~~1~~0\\ 0~~0~~1~~1\\ 0~~0~~0~~1 \end{bmatrix}$，然后算出$P_n,S_n$,则$T_n=nS_n-P_n$，本题结束。

③ P7453 大魔法师

线段树维护矩阵+矩阵乘法。

④ Loj6208 树上询问

线段树+树链剖分+矩阵乘法

⑤ AcWing1305 GT考试

首先, 设$f[i][j]$为到第i个数为止, 最后j个数和不吉利的串t的前缀相同的数量。那么, $f[i-1][j] -> f[i][k]$的话, 就要枚举当前匹配的前缀长度是$j$,第$i$个点取$0-9$的时候, 分别能够转移到哪些$k$

显然, 这个转移是固定的, 可以用$kmp$算出来, 这就说明了$f[i-1]->f[i]$的递推可以有固定的没有变量的转移矩阵，也就是说, 可以用矩阵乘法来做。

假如有$0-9$一共$10$个数, 我们在第$i+1$个位置填$c$, 就可以利用$kmp$求出来此时匹配到了第$j$个位置
也就是说, 当前已经填好了第$i$个数, 准备在第$i+1$的位置填$c$的时候, 可以转移到$j$
$f[i][k] = f[i - 1][0]\times a[0][k] + f[i-1][1]\times a[1][k] + .. + f[m-1][k]\times a[m-1][k]$
我们去累计这样$i->j$的贡献, 就得到了我们的转移矩阵
$X_{i}= X_{i-1} \times A$

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

int n, m, M;
const int N = 20;
int nxt[N + 5];
char s[N + 5];

int mod(int a, int b){
    return (a % b + b) % b;
}

void kmp(){
	for (int i=2, j=0;i<=m;++i){
		while (j && s[j + 1] != s[i]) j = nxt[j];
		if (s[j + 1] == s[i]) ++j;
		nxt[i] = j;
	}
}

void mul(int c[], int a[], int b[][N]){
	int t[N] = {0};
	for (int i=0;i<m;++i){
		for (int j=0;j<m;++j){
			t[i] = (t[i] + a[j] * b[j][i]) % M;
		}
	}
	memcpy(c, t, sizeof(t));
}

void mul(int c[][N], int a[][N], int b[][N]){
	int t[N][N] = {0};
	for (int i=0;i<m;++i){
		for (int j=0;j<m;++j){
			for (int k=0;k<m;++k){
				t[i][j] = (t[i][j] + a[i][k] * b[k][j]) % M;
			}
		}
	}
	memcpy(c, t, sizeof(t));
}

void qmi(int a[][N], int b){
	int t[N][N] = {0};
	for (int i=0;i<m;++i){
		for (int j=0;j<m;++j){
			t[i][j] = (i == j);
		}
	}
	while (b){
		if (b & 1) mul(t, t, a);
		mul(a, a, a);
		b >>= 1;
	}
	memcpy(a, t, sizeof(t));
}

int main(void){
	ios::sync_with_stdio(false);
	cin.tie(0);cout.tie(0);
	cin >> n >> m >> M;
	cin >> s + 1;
	kmp();
	
	int a[N][N] = {0}, f[N] = {1};
	for (int i=0;i<m;++i){
		for (char c='0';c<='9';++c){
			int j = i;
			while (j && s[j + 1] != c) j = nxt[j];
			if (s[j + 1] == c) ++j;
			if (j < m) ++a[i][j];
		}
	}
	qmi(a, n);
	mul(f, f, a);
	int res = 0;
	for (int i=0;i<m;++i){
		res = mod(res + f[i], M);
	}
	cout << res << endl;
	
	return 0;
} 

组合计数

组合数常用计算方法

这两个方法都是嗯……比较抽象和灵活的，在此就以一道例题AcWing1312序列统计的形式简单讲下。

约定：题目中的 $N,L,R$ 用 $n,l,r$ 表示。

寄 有一个写的不错的题解，直接看他的吧，有什么不详细的地方我再补充0.0

1. 映射法

2. 隔板法

隔板法我在基础篇里也有介绍，在这里，隔板法是要求所有盒子非空的，这是应用隔板法的一个前提条件。而在法1里面对于公式的化简，也是一个常见套路：先加上一个常数值，再在最后减去以简化计算。

代码：

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

int n, L, R;
const int M = 1e6 + 3;

int mod(int a, int b){
	return (a % b + b) % b;
}

int qmi(int a, int b, int p){
	int res = 1 % p;
	while (b){
		if (b & 1) res = (ll) res * a % p;
		a = (ll) a * a % p;
		b >>= 1;
	}
	return res;
}

int C(int a, int b, int p){
	if (b > a) return 0;
	int up = 1, down = 1;
	for (int i=1, j=a;i<=b;++i, --j){
		up = (ll)up * j % p;
		down = (ll)down * i % p;
	}
	int res = (ll)up * qmi(down, p - 2, p) % p;
	return res;
}

int lucas(int a, int b, int p){
	if (a < p && b < p) return C(a, b, p);
	return (ll)lucas(a / p, b / p, p) * C(a % p, b % p, p) % p;
}

int main(void){
	int t;
	scanf("%d", &t);
	while (t--){
		scanf("%d%d%d", &n, &L, &R);
		int res = lucas(R - L + n + 1, R - L + 1, M);
		printf("%d\n", mod(res - 1, M));
	}
	return 0;
}

悬线法

基础模板

oi-wiki的定义：https://oi-wiki.org/misc/hoverline/

一个简单的题目，简单理解悬线法的过程：

SP1805 Largest Rectangle in a Histogram

https://www.luogu.com.cn/problem/SP1805

首先，$n$个矩形就相当于$n$条悬线，我们知道最大的面积肯定是由某一个悬线向左右扫过而形成的。那么悬线的扩展，显然的满足某一递推关系，可以帮助我们将复杂度由原来的$O(N^2)$变成$O(N)$。

定义$l_i$为当前的$i$位置能扩展到的悬线的最左端，初始时$l_i=i$。假设已经处理好了前$i-1$个位置的答案，那么当$h_i\leq h_{i-1}$时，$i$也能扩展到$i-1$能扩展到的位置。如果$h_i\leq h_{l_{i-1}-1}$时，又可以接着往前扩展……直到扩展到边界，我们就停止。因此对于$\forall i$，我们有

while (L[i] > 1 && a[i] <= a[L[i] - 1]) L[i] = L[L[i] - 1];

同样的，假如我们已经处理好$i+1$到$n$的答案，$r_i$也能不断的向右扩展。

while (R[i] < n && a[i] <= a[R[i] + 1]) R[i] = R[R[i] + 1];

那么完整代码：

int n;
const int N = 1e5 + 5;
int a[N], L[N], R[N];

void solve(){
	for (int i=1;i<=n;++i){
		scanf("%d", &a[i]);
		L[i] = R[i] = i;
	}
	for (int i=1;i<=n;++i){
		while (L[i] > 1 && a[i] <= a[L[i] - 1]) L[i] = L[L[i] - 1];
	}
	for (int i=n;i>=1;--i){
		while (R[i] < n && a[i] <= a[R[i] + 1]) R[i] = R[R[i] + 1];
	}
	ll res = 0;
	for (int i=1;i<=n;++i){
		res = max(res, 1LL * a[i] * (R[i] - L[i] + 1)); 
	}
	printf("%lld\n", res);
}

UVA1619/POJ2796 Feel Good

给出长度为$n$的数组$a[n]$，找到一个子区间，使得子区间内的最小值与区间内所有元素和的乘积最大，如果有多个答案，输出长度最小的答案，如果仍有多个答案，输出最左端序号最小的答案。

枚举这个最小值，它一旦向左右扩展，就肯定会增加这个乘积的值，这样的话，又变成了一个悬线法求最大子矩形的问题。

数据范围：$1\leq n \leq 10^5$ 我恨UVA的多组数据和格式

int n;
const int N = 1e5 + 5;
int a[N], L[N], R[N];
ll pre[N];
bool fst = 1;

void solve(){
//	n = read();
	if (!fst){
		puts("");
	}
	fst = 0;
	for (int i=1;i<=n;++i){
		a[i] = read();
		pre[i] = pre[i - 1] + a[i];
		L[i] = R[i] = i;
	}
	for (int i=1;i<=n;++i){
		while (L[i] > 1 && a[i] <= a[L[i] - 1]) L[i] = L[L[i] - 1];
	}
	for (int i=n;i>=1;--i){
		while (R[i] < n && a[i] <= a[R[i] + 1]) R[i] = R[R[i] + 1];
	}

	ll res = 0;
	int aL = 1, aR = 1;
	for (int i=1;i<=n;++i){
		ll cur = (pre[R[i]] - pre[L[i] - 1]) * a[i];
		if (cur > res){
			res = cur, aL = L[i], aR = R[i];
		}
		else if (cur == res){
			if (R[i] - L[i] < aR - aL){
				aL = L[i], aR = R[i];
			}
			else if (R[i] - L[i] == aR - aL){
				if (L[i] < aL){
					aL = L[i], aR = R[i];
				}
			}
		}
	}
	printf("%lld\n%d %d\n", res, aL, aR);
}

int main(void){
	while (~scanf("%d", &n)){
		if (n == 0){
			puts("");
		}
		else solve();
	}
	return 0;
}

最大子矩形：p4147 玉蟾宫

嗯……差不多捏，问题。oiwiki留的课后习题也写了，就不放出来了。

---

xcpc真题 King’s Children

有一个$n\times m$的网格，里面有最多$A~to ~Z$个孩子，king要把这个网格分成若干矩形，要满足：

• 每个矩形必须精确的包括一个孩子
• 每个格子都必须精确的属于1个矩形（也就是矩形不相交的分完整个大网格）
• 包括了$A$的矩形的面积尽可能大

数据范围$1\leq n,m\leq 1000$

K题我的思路就是，每个矩形都用悬线法进行选取和填充。但是，由于填充顺序的不同，有极低的概率出现最后有矩形没有被完全填上的情况，因此做一个简单的check，如果填充错误，则随机化顺序，重新填充答案。大概2~3次随机后就不可能出现还没填充的情况了，所以这个复杂度是完全可行的。

代码：

int n, m, sx, sy;
const int N = 1000 + 5;
char s[N][N], ss[N][N];
int U[N], L[N], R[N];
struct node{
	char ch;pii cor;
};
vector<node> alp;

bool check(pii s, pii a, pii b){
	return (s.xx >= a.xx && s.xx <= b.xx) && (s.yy >= a.yy && s.yy <= b.yy);
}

void putin(char ch, pii cor){
	int res = 0; pii r1 = cor, r2 = cor;
	for (int j=1;j<=m;++j) U[j] = 0;
	for (int i=1;i<=n;++i){
		for (int j=1;j<=m;++j){
			if (s[i][j] == '.' || s[i][j] == ch){
				U[j]++;
			}
			else{
				U[j] = 0;
			}
			L[j] = R[j] = j;
		}
		for (int j=1;j<=m;++j){
			while (L[j] > 1 && U[j] <= U[L[j] - 1]) L[j] = L[L[j] - 1];
		}
		for (int j=m;j>=1;--j){
			while (R[j] < m && U[j] <= U[R[j] + 1]) R[j] = R[R[j] + 1];
		}
		for (int j=1;j<=m;++j){
			int cur = U[j] * (R[j] - L[j] + 1);
			pii c1 = pii(i - U[j] + 1, L[j]), c2 = pii(i, R[j]);
			if (cur > res && check(cor, c1, c2)){
				res = cur, r1 = c1, r2 = c2;
			}
		}
	}
	for (int i=r1.xx;i<=r2.xx;++i){
		for (int j=r1.yy;j<=r2.yy;++j){
			if (s[i][j] == '.'){
				s[i][j] = 'a' + (ch - 'A');
			}
		}
	}
}

void solve(){
	for (int i=1;i<=n;++i){
		for (int j=1;j<=m;++j){
			s[i][j] = ss[i][j];
		}
	}
	int alen = alp.size(), t = alen - 1;
	if (alen > 1){
		for (int i=0;i<alen-1;++i){
			int x = 1 + rand() % t;
			swap(alp[x], alp[t]);
			t--;
		}
	}
	
	for (int i=0;i<alp.size();++i){
		putin(alp[i].ch, alp[i].cor);
	}
	
}

bool isOK(){
	for (int i=1;i<=n;++i){
		for (int j=1;j<=m;++j){
			if (s[i][j] == '.'){
				return false;
			}
		}
	}
	return true;
}

int main(void){
	srand(time(NULL));
	int T;
	T = 1;
	ios::sync_with_stdio(false);
	cin.tie(0);cout.tie(0);
	cin >> n >> m;
	for (int i=1;i<=n;++i){
		for (int j=1;j<=m;++j){
			cin >> ss[i][j];
			if (ss[i][j] == 'A'){
				alp.push_back(node{'A', pii(i, j)});
			}
		}
	}
	for (int i=1;i<=n;++i){
		for (int j=1;j<=m;++j){
			if (ss[i][j] != '.' && ss[i][j] != 'A'){
				alp.push_back(node{ss[i][j], pii(i, j)});
			}
		}
	}
	do{
		solve();
	}while(!isOK());
	for (int i=1;i<=n;++i){
		for (int j=1;j<=m;++j){
			cout << s[i][j];
		}
		cout << endl;
	}
	
	return 0;
}

数论分块

简介

数论分块通常被用来以$O(\sqrt n)$的复杂度快速计算形如$\sum \limits_{i=1}^n f(i)g(\lfloor \frac n i \rfloor)$的含有除法向下取整的和式，它的核心思想是将$\lfloor \frac n i \rfloor$相同的数打包同时计算，主要利用了Fubini定理。

证明

1. 证明时间复杂度为$O(\sqrt n)$

**引理1 **对于任意一个正整数$n$，$\lfloor \frac n d \rfloor(d\in[1,n])$的数量级为$\sqrt n$。

• $d\leq \sqrt n$，假设所有的$\lfloor \frac n d \rfloor$取值均不同，则存在$\sqrt n$种结果
• $d>\sqrt n$，则$\lfloor \frac n d \rfloor<\sqrt n$，取值同样有$\sqrt n$。

因此，所有的可能取值一定$\leq 2\sqrt n$，即整数分块的复杂度为$O(\sqrt n)$。

2. 证明算法的正确性

引理2 对于$\forall i,$ 若满足$\lfloor \frac n i \rfloor=C$ ，则所有满足条件的$i$一定是一段连续区间的集合。

反证法，若不是连续区间，则一定$\exist t\in(l,r)$，使得$\lfloor \frac n l \rfloor=\lfloor \frac n r \rfloor$且$\lfloor \frac n l \rfloor \neq \lfloor \frac n t \rfloor$。

$\because t>l,t<r$

$\therefore \lfloor \frac n t \rfloor \geq \lfloor \frac n l \rfloor, \lfloor \frac n t \rfloor \leq \lfloor \frac n r \rfloor$ ，即$\lfloor \frac n t \rfloor=\lfloor \frac n l \rfloor$，矛盾，故得证。

引理3 若满足$\lfloor \frac n i \rfloor=C$的$i$集合为$i\in[l,r]$，则只需要知道$l$，就可以求出对应的$C,r$。

显然，$C=\lfloor \frac n l \rfloor$，下面就要证明$r=\lfloor \frac{n}{\lfloor \frac n l \rfloor} \rfloor=\lfloor \frac n C \rfloor$

令$p=\lfloor \frac n C \rfloor$，则$n=p\times C+k,k\in[0,\min\{p,C\})$，下面要证明$p$是满足$i$的性质的最大的整数。

• $p$满足$i$的性质，即$\lfloor \frac n p \rfloor=C$

$p=\frac{n-k}{C}$，则$\frac n p=\frac{C\times n}{n-k}\geq \frac {C\times n}{n}$，即$\frac n p \geq C$。

又$n\mod p=k,k<p$，则$n-k<p\times (C+1)$，即$\lfloor \frac n p \rfloor<C+1$

因此$\lfloor \frac n p \rfloor=C$

• $p$是使得$\lfloor \frac n i \rfloor=C$成立的最大的数。

反证法，假如$p$不是最大的数，则$p+1$一定可以使得$\lfloor \frac n {p+1} \rfloor=C$

$\therefore (p+1)\times C=n-g(g\in[0,\min\{p+1,C\})$

$\therefore p\times C + C=n-g$，又$p\times C=n-k$

所以$C=k-g$，又$k,g\in[0,C)$，所以矛盾，得证。

综上，$\lfloor \frac n C \rfloor$即为令$\lfloor \frac n i \rfloor=C$的最大的数

模板

ll H(ll n){
	ll res = 0, l = 1, r;
	while (l <= n){
		r = n / (n / l);
		res = res + (r - l + 1) * (n / l);
		l = r + 1;
	}
	return res;
}

例题

P2261余数之和

给出$n,k$，计算$G(n,k)=\sum\limits _{i=1}^n k \bmod i$

数据范围：$1\leq n,k\leq 10^9$

（懒得自己写题解了，直接cpy洛谷题解orz）由题意得：$ans=∑_{i=1}^nk\bmod i$

我们知道，$a \bmod b = a - b \times \lfloor \frac{a}{b} \rfloor$

因此，$ans = \sum_{i=1}^{n} k - i \times \lfloor \frac{k}{i} \rfloor = nk - \sum_{i=1}^{n} i \times \lfloor \frac{k}{i} \rfloor$

首先枚举块的左边界 $l$，并根据左边界和$ k$ 计算出右边界 $r$。

令 $t = \lfloor \frac{k}{l} \rfloor$，分两种情况讨论：

• $t \neq 0$，则 $r = \min (\lfloor \frac{k}{t} \rfloor , n)$；
• $t = 0$，则 $r = n$。

右边界有了，每一块的和也就可以计算出了。

每一块的和 $=$ 当前块的 $t\times$ 当前块元素个数 $\times$ 当前块 $i$ 的平均值 $= t \times (r-l+1) \times (l+r) \div 2$

当前块处理完后，令 $l = r + 1$，开始计算下一块，直到计算至 n。

ll n, m;

void solve(){
	n = read(), m = read();
	ll res = n * m, l = 1, r;
	while (l <= n){
		if (m / l == 0) r = n;
		else r = min(n, m / (m / l));
		res = res - (r - l + 1) * (m / l) * (l + r) / 2;
		l = r + 1;
	}
	printf("%lld\n", res);
}

数论分块往往和莫比乌斯反演结合起来考察，因此进阶的部分放在可能会存在的莫反专题中一起讲解。

常见数论定义

数论函数

**定义：**一个定义在整数集合上的实数或复数函数$f(n)$，称为一个数论函数，又叫算术函数。：

**举例：$pot _p (n),\mu(n),\varphi(n),\pi(n)$**都是数论函数，$n(i)=n$：一直返回同一个常数$n$。

加性函数与积性函数

定义

• 对于一个数论函数$f:N\rightarrow C$，若$\forall m,n$，满足$gcd(m,n)=1$，有$f(m\times n)=f(m)+f(n)$

则称$f$为加性的，这时定义给出$f(1)=0$。

• 若$\forall m,n$，满足$gcd(m,n)=1$，有$f(m\times n)=f(m)\times f(n)$则称$f$为积性的，这时定义给出$f(1)=1$。

​ 若去掉$gcd(m,n)=1$这一条件仍然满足，则称$f$为完全积性函数。

判定函数为积性函数或完全积性函数

设$n$的唯一分解为$n= \prod_{i=1}^kp_i^{a_i} = p_1^{a_1}p_2^{a_2}..p_k^{a_k}$，那么

• $f(1)=1$且$f(n)=\prod_{i=1}^k f(p_i^{a^i}) \Leftrightarrow f(n)$是积性函数（充要条件）

• $f(1)=1$且$f^{a_1}(p_1)\times f^{a_2}(p_2)\times..\times f^{a_k}(p_k) \Leftrightarrow f(n)$是完全积性函数（充要条件）

这表明，一个积性函数完全由它在素数幂$p_i^{a_i}$上的取值所确定；而完全积性函数则完全由它在素数$p_i$上的取值所确定，我们由此可构造积性函数。

举例：常见的积性函数

单位元函数：$e(n)=[n=1]$，即$e(n)=\begin{cases} 1 & n=1\\ 0 & n\neq 1 \end{cases}$

它卷积上任意的数论函数仍然为原来的数论函数，即满足$f\times e=e\times f=f$

幂函数：$id^k(n)=n^k$，完全积性

单位函数：$id(n)=n$，完全积性，相当于$id^1$

恒等函数：$I(n)=1$（也就是常数值函数，只是刚好这个常数等于1），相当于$id^0$

欧拉函数：$\varphi(n)=\sum\limits _{i=1}^n[gcd(n,i)=1]\times 1$

除数函数：$\sigma_k(n)=\sum\limits _{d|n}d^k$，表示$n$的约数的$k$次幂和

约数和函数：$\sigma(n)=\sigma_1(n)=\sum\limits _{d|n}d$，表示$n$的约数之和

约数个数函数：$\tau(n)=\sigma_0(n)=\sum\limits _{d|n}1$，表示$n$的约数的个数，一般也写作$d(n)$

莫比乌斯函数：$\mu(n)$

生成函数

狄利克雷卷积

定义

1. 狄利克雷卷积(Dirichlet Product)

设现在有两个数论函数$f(n),g(n)$，那么它们的狄利克雷卷积（也叫狄利克雷乘积）也是一个数论函数。记它们的狄利克雷卷积为$h(n)$，则有：

$h(n)=\sum\limits _{d|n}f(d)g(\frac n d)$，或者说$h(n)=\sum\limits _{xy=n}f(x)g(y)$

简记为$h(n)=f(n) * g(n)$

2. 狄利克雷逆(Dirichlet inverse)

若$f*g=e$，则称$g$是$f$的狄利克雷逆，记作$f^{-1}$。

计算狄利克雷逆

首先，当$n=1$时

$(f*f^{-1})(1)=\sum\limits _{d|1}f(d)f^{-1}(\frac 1 d)=f(1)f^{-1}(1)$

又$(f*f^{-1})(1)=e(1)=1$，所以$f(1)f^{-1}(1)=1$

$f^{-1}(1)=\frac 1 {f(1)}$

这说明了，$f^{-1}$存在的必要条件是$f(1)\neq 0$。

当$n>1$时，有

$(f*f^{-1})(n)=\sum\limits _{d|n}f(d)f^{-1}(\frac n d)=f(1)f^{-1}(n)+\sum\limits _{d|n,d>1}f(d)f^{-1}(\frac n d)$

$\because (f*f^{-1})(n)=e(n)=0$

$\therefore f(1)f^{-1}(n)+\sum\limits _{d|n,d>1}f(d)f^{-1}(\frac n d)=0$

$f^{-1}(n)=-\frac1 {f(1)}\sum\limits _{d|n,d>1}f(d)f^{-1}(\frac n d)$

综上，即可递归计算狄利克雷逆

$f^{-1}(n)=\begin{cases} \frac 1 {f(1)} & n=1 \\ -\frac1 {f(1)}\sum\limits _{d|n,d>1}f(d)f^{-1}(\frac n d) & n > 1 \end{cases}$

性质

交换律

$(f * g)(n)=\sum \limits_{xy=n}f(x)g(y)=\sum\limits _{yx=n}=g(y)f(x)=(g*f)(n)$

结合律

分配律

若$f,g$为积性函数，则$h=f*g$也是积性函数

若$g, h=f*g$是积性函数，则$f$也是积性函数

数论函数的卷积关系

0. 对于任意数论函数$f$和恒等函数$I$，有

$(f*I)(n)=\sum\limits_{d|n} f(d)I(\frac n d)=\sum\limits _{d|n}f(d)$

1. 幂函数与除数函数

$(id^k*I)(n)=\sum\limits _{d|n}id^k(n)=\sum\limits n^k=\sigma_k(n)$

即$id^k * I=\sigma_k$。

2. 欧拉函数与恒等函数

$(\varphi*I)=\sum\limits_{d|n}\varphi(d)=n=id(n)$

3. 恒等函数与恒等函数

$I *I=d$（在莫比乌斯函数里证明过了）

莫比乌斯反演

莫比乌斯函数

1. 定义

由唯一分解定理，可以将正整数$n$写成$n= \prod_{i=1}^kp_i^{a_i} = p_1^{a_1}p_2^{a_2}..p_k^{a_k}$的形式，莫比乌斯函数$\mu(n)$的定义为

$\mu(n)=\begin{cases} 1 & n=1 \\ 0 & \exist i, a_i\geq 2\\ (-1)^{k} & \forall i,a_i=1 \end{cases}$

2. 性质

性质1

$\sum\limits _{d|n}\mu(d)=\begin{cases} 1 & n=1 \\ 0 & n \neq 1 \end{cases}$

证明：设$d$为$n$的约数，则$d=\prod_{i=1}^kp_i^{b_i}$，其中$0\leq b_i\leq a_i$。

对于$\mu(d)$，如果$\exist b_i\geq 2$，则$\mu(d)=0$。因此，有贡献的$\mu(d)$一定为$C_k^i\times(-1)^i$，也就是每个质数最多取一次。

则$\sum\limits _{d|n}\mu(d)=\sum\limits _{i=0}^kC_k^i\times(-1)^i$，又$(a-b)^k=\sum\limits_{i=0}^k C_k^i a^kb^{k-i}$

$(1-1)^k=\sum\limits _{i=0}^kC_k^i\times (-1)^k$，故$\sum\limits _{d|n}\mu(d)=0^k=0$

3. 与其他数论函数的关系

(1) $\mu * I = e$

证明：设$n=\prod_{i=1}^kp_i^{a_i}, n'=\prod_{i=1}^kp_i$

则$(\mu*I)(n)=\sum\limits _{d|n}\mu(d)=\sum\limits _{d|n'}\mu(d)\\= \sum\limits _{i=0}^k(-1)^i$

呃，等等，好像性质一已经证明过了啊。$(\mu*I)(n)=[n=1]=e$，

因此，$\mu$是$I$的狄利克雷逆。

(2) $\mu * id = \varphi$

这个在基础篇的性质证明过了QWQ，不写辣

(3) $\mu * d=I$

证明：$(I*I)(n)=\sum\limits _{d|n}I(d)=\sum\limits _{d|n}1=d(n)$

$\therefore d=I*I$，又$\mu=I^{-1}$

$\therefore \mu * d=I$

4. 线性筛法求莫比乌斯函数

void Mobius(int n){
	mu[1] = 1;
	for (int i=2;i<=n;++i){
		if (!st[i]) p[++cnt] = i, mu[i] = -1;
		for (int j=1;p[j]<=n/i;++j){
			st[p[j] * i] = true;
			if (i % p[j] == 0) break;
			mu[p[j] * i] = -mu[i];
		}
	}
}
// 当i为质数时, 显然mu[i]=-1
// 当p[j]为i的最小质数时, 就说明p[j]这个质数出现了>1次, 因此mu[i * p[j]] = 0
// 否则
// (1) mu[i]=0, mu[p[j] * i] = 0
// (2) mu[i]不为0, p[j] * i就相当于增加了一个质数, 因此mu[p[j] * i] = -mu[i]

莫比乌斯反演

莫反的函数定义和转换过程大多依靠平时积累，见过类似套路，就会，没见过，就寄。——yxc

1. 定义

设$f(n)$为数论函数(定义在正整数集合上的函数)

因数形式：

• $F(n)=f*I=\sum\limits_{d|n}f(d) \Leftrightarrow f(n)=\sum\limits _{d|n}\mu(d)\times F(\frac n d)$，

证明（利用狄利克雷卷积）：因为$F(n)=f*I$，则$f=F*I^{-1}=F*\mu$

即$f(n)=\sum\limits_{d|n}\mu(d)\times F(\frac n d)$。

证明（利用性质1+二重积分交换次序的思想）：

$\sum\limits_{d|n}\mu(d)\times F(\frac n d)=\sum\limits_{d|n}\mu(d)\times \sum\limits_{i|\frac n d}f(i)=\sum\limits_{i|n}f(i)\sum\limits_{d|\frac n i}\mu(d)$

（$i$能取到所有$d$可以取到的取值，这样反过来看，把$i$提到前面）

又当且仅当$n=i$时，$\sum\limits_{d|\frac{n}i}\mu(d)=1$，因此$\sum\limits_{d|n}\mu(d)\times F(\frac n d)=f(n)$

倍数形式：

• $F(n)=\sum\limits_ {n|N}f(N) \Leftrightarrow f(n)=\sum\limits_{n|N}F(N)\mu(\frac N n)$，（枚举$N$为$n$的所有倍数，$N\in[n,+\infin)$）

证明：$\sum\limits_{n|N}F(N)\mu(\frac N n)=\sum\limits_{n|N}\mu(\frac N n)\sum\limits _{N|i}f(i)$

设$d=\frac N n$，则$N=dn$，则$dn|i$，即$d|\frac i n$

因此$\sum\limits_{n|N}\mu(\frac N n)\sum\limits _{N|i}f(i)=\sum\limits_{d|\frac i n}\mu(d)\sum\limits _{N|i}f(i)$

又当且仅当$n=i$时，$\sum\limits_{d|\frac{n}i}\mu(d)=1$，因此$f(n)=\sum\limits_{n|N}F(N)\mu(\frac N n)$

运用莫反的时候，通常都是因为$F(n)$好求，但是$f(n)$不好求，因此将$f(n)$用$F,\mu$表示出来。

2. 应用1：莫反+整数分块

p2522 Problem b

数据范围：$1\leq n,k\leq 5\times 10^4;1\leq a\leq b\leq 5\times 10^4;1\leq c \leq d \leq 5\times 10^4$

思路：详细的整理一下吧。

首先，题目要我们求的东西，可以先拆成一个二维前缀和，$A[a,b][c,d]=A[1,b][1,d]-A[1,b][1,c-1]-A[1,a-1][1,d]+A[1,a-1][1,c-1]$。

设$f(k)=\sum\limits _{x=1}^a\sum\limits _{y=1}^b[(x,y)=k]$，然后我们方便求的是这个$F(k)=\sum\limits _{x=1}^a\sum\limits _{y=1}^b[k|(x,y)]$，且$F(k)=\sum\limits _{k|N}f(N)$

则代入莫反倍数形式得$f(k)=\sum\limits _{k|N}\mu(\frac N k) F(N)$

先求$F(N)$。首先，$N|(x,y)$，也就是说，$N|x,N|y$，因此所有满足条件的点对数量为$\lfloor \frac a N \rfloor\times \lfloor \frac b N \rfloor$

则$f(k)=\sum\limits _{k|N}\mu(\frac N k)\lfloor \frac a k \rfloor\times \lfloor \frac b k \rfloor$，设 $t=\frac{N}{k}$，显然枚举 $t$的结果为$1,2,..,$这样的整数，$N=tk$。

$f(k)=\sum\limits_{t}\mu(t)\lfloor \frac a {tk} \rfloor\times \lfloor \frac b {tk} \rfloor$，再运用整数分块的知识进行求解即可，注释都写在代码里吧。

#include <bits/stdc++.h>
using namespace std;
#define ll long long
typedef pair<int, int> pii;
typedef pair<ll,ll> pll;
#define xx first
#define yy second
#define ls (oo << 1)
#define rs (oo << 1 | 1)
#define PI acos(-1.0)

ll read(void);

int n, cnt;
const int N = 5e4 + 5; 
int p[N], mu[N];
int pre[N];
bool st[N];

//求Mobius函数和前缀和(分块的时候用)
void Mobius(int n){
	mu[1] = 1;
	for (int i=2;i<=n;++i){
		if (!st[i]) p[++cnt] = i, mu[i] = -1;
		for (int j=1;p[j]<=n/i;++j){
			st[p[j] * i] = true;
			if (i % p[j] == 0) break;
			mu[p[j] * i] = -mu[i];
		}
	}
	for (int i=1;i<=n;++i){
		pre[i] = pre[i - 1] + mu[i];
	}
}

ll f(int a, int b, int k){
    a /= k, b /= k;
	ll res = 0, n = min(a, b), l = 1, r;
    // 在[l,r]这段,(a/l)*(b/l)为定值,那么展开和式, 可以打包计算这一部分的和为(定值*mu的前缀和)
	while (l <= n){
		r = min(n, min(a / (a / l), b / (b / l)));
		res += 1LL * (pre[r] - pre[l - 1]) * (a / l) * (b / l);
		l = r + 1;
	}
	return res;
}

void solve(){
	int a, b, c, d, k;
	a = read(), b = read(), c = read(), d = read(), k = read();
    // 二维前缀和,或者说一个简单的容斥
	ll res = f(b, d, k) - f(b, c - 1, k) - f(a - 1, d, k) + f(a - 1, c - 1, k);
	printf("%lld\n", res);
}

int main(void){
	int T;
	Mobius(N - 1);
	T = read();
	while (T--){
		solve();
	}
	
	return 0;
}

ll read(void){
    ll x = 0, f=1;char ch;
    do{ch = getchar();if (ch == '-') f=-1;}while(ch<'0' || ch>'9');
    do{x = x*10 + (ch-'0');ch = getchar();}while(ch>='0' && ch<='9');
    return x*f;
}

/*
敬告kz: 
====================================
  1. 相信自己 
  2. 看清题意, 考虑清楚再动手 
  3.   **** 今天的数组有没有开小呀 ？ ****  **** 今天的数组有没有开小呀 ？ ****
  4. 是不是想复杂了？ 
  5. 数据溢出？
  6. 数组越界？边界情况？ 
  6. 不要犯低级错误！！！ 时间复杂度？空间复杂度？精度有没有问题？ 
====================================
* 提交的时候注意看编译器！c++17 / c++20 / python3 
*/ 

3. 应用2：莫反+提取公因数

p3327约数个数和 莫反+双分块

设$d(x)$为$x$的约数个数，给定$T$组$n,m$，求$\sum\limits _{i=1}^N \sum\limits_{j=1}^M d(i\times j)$

数据范围：$1\leq N,M,T\leq 5\times 10^4$

$\sum\limits _{i=1}^N \sum\limits_{j=1}^M d(i\times j)=\sum\limits _{i=1}^N \sum\limits_{j=1}^M \sum\limits _{x|i} \sum\limits_{y|j} [(x,y)=1]$

证明：设$i=\prod_{i=1}^k p_i^{a_i},j=\prod_{i=1}^k p_i^{b_i}$，$0\leq a_i,b_i$

则$i\times j=\prod_{i=1}^k p_i^{a_i+b_i}$，$d(i\times j)=\prod_{i=1}^k(a_i+b_i+1)$

即从$i$中选出约数$x$，$j$中选出约数$y$，对于$p_1$而言，若要求$(x,y)=1$

则可以$x=1,y=1$，或者$x=1,y=\in[p_1,p_1^{b_1}]$，或者$x\in[p_1,p_1^{a_1}],y=1$

一共是$(a_1+b_1+1)$种取法，其他质数同理。根据乘法原理，这些取法正好就是$d(i\times j)$。

• 设出$f(n),F(n)$。

设$f(n)=\sum\limits _{i=1}^N \sum\limits_{j=1}^M \sum\limits _{x|i} \sum\limits_{y|j} [(x,y)=n]$，显然$f(1)$就是答案。

设$F(n)=\sum\limits _{i=1}^N \sum\limits_{j=1}^M \sum\limits _{x|i} \sum\limits_{y|j} [n|(x,y)]$，则$F(n)=\sum\limits _{n|d}f(d)$

即$f(n)=\sum\limits _{n|d}\mu(\frac d n)F(d)$令$T=\min(N,M)$，则$f(1)=\sum\limits _{d=1}^T\mu(d)F(d)$。

• 再化简$F$。

$F(n)=\sum\limits _{i=1}^N \sum\limits_{j=1}^M \sum\limits _{x|i} \sum\limits_{y|j} [n|(x,y)]=\sum\limits _{x=1}^N \sum\limits_{y=1}^M \lfloor \frac N x \rfloor \lfloor \frac M y \rfloor [n|(x,y)]$

证明：首先，$x|i,y|j$，那么$x,y$肯定是能取到$[1,N],[1,M]$的。当$x,y$固定后，$[n|(x,y)]$和$i,j$是没有关系的，我们可以把它提出来。那么，里面就变成了$\sum\limits _{i=1}^{\lfloor \frac N x \rfloor}\sum\limits _{j=1}^{\lfloor \frac M y \rfloor}1$，也就是$N,M$里面有多少个$i,j$，它们是$x,y$的倍数，得证。

下面再消掉$[n|(x,y)]$这个条件。

设$x'=\lfloor \frac x n \rfloor,y'=\lfloor \frac y n \rfloor$

$F(n)=\sum\limits _{x=1}^N \sum\limits_{y=1}^M \lfloor \frac N x \rfloor \lfloor \frac M y \rfloor [n|(x,y)]=\sum\limits _{x'=1}^{\lfloor \frac N n \rfloor}\sum\limits _{y'=1}^{\lfloor \frac M n \rfloor}\lfloor \frac N {nx'} \rfloor\lfloor \frac M {ny'} \rfloor$

令$N'=\lfloor \frac N n \rfloor,M'=\lfloor \frac M n \rfloor$

$F(n)=\sum\limits _{x'=1}^{N'} \sum\limits_{y'=1}^{M'} \lfloor \frac {N'} {x'} \rfloor \lfloor \frac {M'} {y'} \rfloor=(\sum\limits _{x'=1}^{N'} \lfloor \frac {N'} {x'} \rfloor)\times(\sum\limits_{y'=1}^{M'} \lfloor \frac {M'} {y'} \rfloor)$

令$h(n)=\sum\limits_{i=1}^{n} \lfloor \frac {n} {i} \rfloor)$，也就是标准整数分块，则$F(n)=h(N')\times h(M')$。

• 再求$f(1)$

$f(1)=\sum\limits _{d=1}^T\mu(d)h(\lfloor \frac N d \rfloor)h(\lfloor \frac M d \rfloor)$

由于$h(x)$只和$x$有关，所以可以再分一次块，因此每次查询复杂度$O(\sqrt N)$，总时间复杂度$O(N\sqrt N)$。

int cnt;
const int N = 5e4 + 5;
int p[N], h[N], pre[N], mu[N];
bool st[N];

void Mobius(int n){
	mu[1] = 1;
	for (int i=2;i<=n;++i){
		if (!st[i]) p[++cnt] = i, mu[i] = -1;
		for (int j=1;p[j]<=n/i;++j){
			st[p[j] * i] = true;
			if (i % p[j] == 0) break;
			mu[p[j] * i] = -mu[i];
		}
	}
	for (int i=1;i<=n;++i){
		pre[i] = pre[i - 1] + mu[i];
	}
}

void H(int n){
	for (int i=1;i<=n;++i){
		for (int l=1, r;l<=i;l=r + 1){
			r = min(i, i / (i / l));
			h[i] += (r - l + 1) * (i / l); 
		}
	}
}

void solve(){
	int n, m;
	n = read(), m = read();
	ll res = 0;
	int k = min(n, m);
	for (int l=1, r;l<=k;l=r + 1){
		r = min(k, min(n / (n / l), m / (m / l)));
		res += (ll)(pre[r] - pre[l - 1]) * h[n / l] * h[m / l];
	}
	printf("%lld\n", res);
}

int main(void){
	int T;
	Mobius(N - 1);
	H(N - 1);
	T = read();
	while (T--){
		solve();
	}
	
	return 0;
}

丑数筛

给定一个质数集合$S=\{p_1,p_2,..,p_k\}$，只由这些质数相乘得到的数我们成为丑数（Humble/Ugly Numbers）。习惯上，我们认为第一个丑数是1，但是也可能不是，所以看清题意。记第$n$大 的丑数为$h[n]$。

无论哪种方法，这个丑数实际上都是可能非常大的，建议直接上__int128！！

第一种筛法，也是最常用的，就是搞一个优先队列，每次提出来优先队列中最小的数，注意去重。

这种方法的时间复杂度是$O((n\times k)\log (n\times k))$

#define ll __int128
const ll INF = 9e37, N = 5e4 + 5;
ll h[N], p[3] = {2, 3, 5}, cnt[3];
priority_queue<ll, vector<ll>, greater<ll> > pq;
void Humble(int n){
	pq.push(1);
    ll last  = 1;
    for (int i=1;i<=n;++i){
        for (int j=0;j<3;++j){
            pq.push(last * p[j]);
        }
        last = pq.top();
        while (!pq.empty() && last == pq.top()) pq.pop();
        h[i] = last;
    }
}

第二种线性筛法，相当于做了一个DP+双指针，复杂度$O(n\times k)$。注意，这里的下标是从0开始的，并且，INF要足够足够大！！！！

ll h[N], p[3] = {2, 3, 5}, cnt[3];
void Humble(int n){
	h[0] = 1; // 下标0的地方是第一个数
	for (int i=1;i<=n;++i){
		h[i] = INF;
		for (int j=0;j<3;++j){
			while (p[j] * h[cnt[j]] <= h[i - 1]) cnt[j]++;
			if (p[j] * h[cnt[j]] < h[i]){
				h[i] = p[j] * h[cnt[j]];
			}
		}
	}
}
// cnt[j]记录的是当前乘的最后一个质因子是p[j]的, 最小的那个数的下标
// h[i]一定从这些数中产生

快速幂&快速乘

快速幂

1. 正常的快速幂

// 求a^b % p 时间复杂度O(log b)
ll qmi(ll a, ll b, ll p){
    ll res = 1LL % p;
    while (b){
        if (b & 1) res = res * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return res;
}

2. 欧拉降幂

运用基础篇欧拉定理推论1公式：$A^b\equiv A^{b\mod \varphi(m)} (\mod m)$，即可。因为经常用到高精度，所以在这里放一个Python的代码吧

# -*- coding: utf-8 -*-
"""
Created on Mon Apr 11 14:47:12 2022
@ 求x^(n-1) % mod
@author: KZ
"""

def Eulers(x):
    res = x
    for i in range(2, x):
        if i > x // i:
            break
        if x % i == 0:
            res = res // i * (i - 1)
            while x % i == 0:
                x = x // i
    if x > 1:
        res = res // x * (x - 1)
    return res

def qmi(a, b, p):
    res = 1 % p
    while (b):
        if (b & 1 == 1):
            res = res * a % p
        a = a * a % p
        b = b >> 1
        
    return res

arr = input().split()
x = int(arr[0])
n = int(arr[1])
mod = int(arr[2])
if mod == 1:
    print(0)
elif n == 1:
    print(1)
else:
    M = Eulers(mod)
    b = (n - 1) % M
    res = qmi(x, b, mod)
    print(res)

3. 同一底数和同一模数的预处理快速幂

前置知识：数论分块；数学基础篇相关定理

4. 十进制快速幂

当二进制快速幂都超时的时候，可以尝试用下十进制快速幂。

比如，$3^{405}=(3^1)^5\times(3^{10})^0+(3^{100})^4$

感觉不太能用到的东西（欧拉降幂和python高精都可以解决吧应该），随便贴个代码参考

#include<iostream>
#include<cstdio>
#include<cstring>
using namespace std;
#define ri register int 
#define ll long long 
ll b,t,h,P;

char c[ 100007 ];
ll TenthPow( ll a ){
    ll ans = 1,s = a;
    while( t >= 0 ){
        ll cnt = c[ t ] - '0',cur = s;
        for( ri i = 1 ; i <= cnt ; ++i )
          ans = ans * s % P;
         //和二进制快速幂不同的地方之一，请结合上面列举的拆分过程理解
        for( ri i = 1 ; i < 10 ; ++i )
          cur = cur * s % P;
        //进位
        s = cur;
        ans %= P;
        --t;
    }
    return ans;
}
int main(){
	cin >> b >> c >> P;
	t = strlen( c );
    --t;//字符串读入指数，指数可能达到几十万位
	cout << b << "^" << c << " mod " << P << "=" << TenthPow( b );
	return 0;
}

快速乘

1. 虚假的快速乘 复杂度$O(\log b)$

// 求a^b % p 时间复杂度O(log b)
ll qmul(ll a, ll b, ll p){
    ll res = 1LL % p;
    while (b){
        if (b & 1) res = (res + a) % p;
        a = (a + a) % p;
        b >>= 1;
    }
    return res;
}

2. 真正的快速乘 复杂度$O(1)$

应用范围：计算$a\times b \bmod m;a,b\leq m\leq 10^{18}$

等，等一下，$a,b\leq 10^{18}$，但是$\_\_int128$极限可以处理$10^{37}$这么大的数，为什么还要快速乘？只有在数组空间开不下的很苛刻的情况下才能用得上吧。那这里就截图贴一下推导过程和代码。

#define ull unsigned long long
#define ll long long
#define ld long double
ll binmul(ll a, ll b, ll m) {
  ull c =
      (ull)a * b -
      (ull)((ld)a / m * b + 0.5L) * m;
  if (c < m) return c;
  return c + m;
}

不得不说的那些题

拆数+数学思维 102780J Something that resembles Waring’s problem

给定$N$，要求构造出$x_1^3+x_2^3+...+x_k^3=N$，$k\leq 5$

数据范围：$1\leq N\leq 10^{100000};|x_i|\leq 10^{110000}$

这道题，应当通过立方和之间的加减，消去3次和2次项，尝试构造出一次项。通过神秘的东方力量日常的积累，也许可以蒙出来想到$(x-1)^3+(x+1)^3-x^3-x^3=6x$

认真的说，这个更像完全平方公式的拓展：$(x+1)^2-(x-1)^2-x^2-x^2=4x$

总之，我们可以通过这个公式，构造出一个$6$的倍数，然后让这个数加某一个$t^3$，使得它变成$N$。

也就是说，$N=6x+t^3$。进一步的，我们发现这样一个神奇性质：

$1^3\bmod 6=1;2^3\bmod 6=2;3^3\bmod 6=3;4^3\bmod 6=4;5^3\bmod 6=5;$

因此我们发现，令$t=N \bmod 6$，让$a=\lfloor \frac N 6\rfloor-\lfloor \frac t 6\rfloor$

则$N=6\times \lfloor \frac N 6\rfloor+N\bmod 6=6x+t\bmod 6=6x-6\lfloor \frac t 6\rfloor+t^3=6a+t^3$。

Python代码如下：

n = int(input())
t = n % 6
a = (n // 6) - (t**3) // 6
print(5)
print("{0} {1} {2} {3} {4}".format(a-1, a+1, -a, -a, t))

数学+单调栈 102483A Access Points

给定$n$个点，可以改变它们的坐标位置，对于每一个点$(x,y)$，将$(x,y)$改变为$(x',y')$的花费为$(x'-x)^2+(y'-y)^2$。经过若干次改变之后，它们必须要满足$\forall i<j$，点$i,j$的坐标分别为$(x_i,y_i),(x_j,y_j)$，$x_i\leq x_j,y_i\leq y_j$。问最小的总花费之和为多少。

数据范围：$1\leq n\leq 10^5$，每个点坐标$(s_i,t_i)$，$1\leq s_i,t_i\leq 10^6$。

首先，$x,y$方向的花费其实是相互独立的，那么答案就是$\sum\limits _{i=1} ^n (x'-x)^2$与$\sum\limits _{i=1} ^n (y'-y)^2$之和，独立计算。

问题转化为，对于一个数组$a[n]$，可以对所有的$a_i$执行加或减操作，变成$p_i$，求让$p[n]$单调非递减的最小花费。首先，如果$a[n]$本身就是单调非递减的，那么不需要进行任何操作，$\forall i\in[1,n]$，$p[i]=a[i]$，总花费为0。如果$a[n]$不满足该规律，可能有增有减，那么可以考虑把$a[n]$分割为若干个串，对于每个串，都将这一个串中的数变成它们的均值，再处理下一个串，使得下一个串的平均值大于上一个串的平均值。也就是每个串的均值递增，串内部的所有数都变为相等。

首先，这样显然是可以构造出一个非递减序列的。然后我们来证明它可以得到最小花费。

如果某个$a_i\leq a_{i+1}$，它是不需要操作的，只有当$a_i> a_{i+1}$的时候才需要操作。那么，让$p[i]=p[i+1]$的花费，一定比让$p[i]<p[i+1]$的花费要更小。（邻项微扰法就可以反证出）那么对于某个值，它要么会和前面若干个数合并，要么会和后面若干个数合并（这个值本身就合法时也可以保持自己一人一组）。如果它的值比前面的$p[i-1]$所属组的均值要小或者相等，那么它应该归到$p[i-1]$组，因为规定了$p[i+1]$组的均值大于$p[i-1]$。否则，应该归到$p[i+1]$组（因为这不会带来花费）。用一个单调栈就可以很方便的处理。

and一组数变为一个值的时候，为什么是变成均值最小呢？这是初中数学啊QAQ

让我们回顾一下最小二乘法$\epsilon=\sum(y-y_i)^2$，$\epsilon$最小时，取到真值。对这个二次函数求导

$\frac{d}{dy}\epsilon=\frac{d}{dy}\sum(y-y_i)^2=2\sum(y-y_i)=2(n\times y-\sum y_i)=0$

显然当$y=\frac{\sum y_i}{n}$时取到最小值。

int n, tt;
const int N = 1e5 + 10;
ll x[N], y[N];
pll s[N];
double ans;

double cal(ll a[]){
	tt = 0, ans = 0;
	for (int i=1;i<=n;++i){
		ll cur = a[i], le = 1;
		while (tt && s[tt].xx*le>=s[tt].yy*cur){
			le += s[tt].yy, cur += s[tt].xx;
			--tt;
		}
		s[++tt] = pll(cur, le);
	}
	int inx = n;
	while (tt){
		double tmp = 1.0*s[tt].xx/s[tt].yy;
		ll le = s[tt].yy; --tt;
		for (int j=0;j<le;++j){
			ans += (tmp-a[inx-j]) * (tmp-a[inx-j]);
		}
		inx -= le;
	}
	return ans;
}

int main(void){
	n = read();
	for (int i=1;i<=n;++i){
		x[i] = read(), y[i] = read();
	}
	printf("%.9lf\n", cal(x)+cal(y));
}

概率DP+状压 103447G Damaged Bicycle

你需要从1 走到n ，初始速度是t，某些地方有自行车，每个位置自行车有pi的概率是坏掉的，如果自行车没坏可以骑上自行车，速度是r，可以一直骑着到终点。

懒得写题解了0.0，记录下Dijstra+状压DP模板

int n, m, t, r, k;
const int N = 1e5 + 5, M = 20 + 5, K = 1e6 + 5;
const int INF = 2e9 + 10;
int dis[M][N], cnt[N], b[M], id[N];
bool vis[M][N];
double f[M][K], p[M];
vector<pii> ve[N];
priority_queue<pii, vector<pii>, greater<pii> > pq;

void Dijkstra(int id, int u){
	for (int i=1;i<=n;++i) dis[id][i] = INF;
	pq.push(pii(0, u)); dis[id][u] = 0;
	while (!pq.empty()){
		pii tmp = pq.top(); pq.pop();
		int v = tmp.yy, w = tmp.xx;
		if (vis[id][v]) continue;
		vis[id][v] = true;
		dis[id][v] = w;
		for (auto i : ve[v]){
			if (dis[id][i.yy] > dis[id][v] + i.xx){
				dis[id][i.yy] = dis[id][v] + i.xx;
				pq.push(pii(dis[id][i.yy], i.yy));
			}
		}
	}
}

double dp(int st, int u){
	if (f[u][st] > 0) return f[u][st];
	double time = p[u] * dis[u][n] / t + (1.0 - p[u]) * dis[u][n] / r;
	for (int i=0;i<k;++i){
		if (st >> i & 1) continue;
		time = min(time, p[u] * (1.0 * dis[u][b[i + 1]] / t + dp(st | (1<<i), i+1)) + (1.0-p[u]) * dis[u][n] / r);
	}
	return f[u][st] = time;
}

int main(void){
	scanf("%d%d", &t, &r);
	scanf("%d%d", &n, &m);
	int u, v, w, cp;
	for (int i=1;i<=m;++i){
		scanf("%d%d%d",&u, &v, &w);
		ve[u].push_back(pii(w, v));ve[v].push_back(pii(w, u));
	}
	scanf("%d", &k);
	for (int i=1;i<=k;++i){
		scanf("%d%d", &b[i], &cp);
		p[i] = cp / 100.0, id[b[i]] = i;
		Dijkstra(i, b[i]);
	}
	Dijkstra(k + 1, 1), Dijkstra(k + 2, n);
	if (dis[k+1][n] >= INF){
		puts("-1"); return 0;
	}
	p[k + 1] = 1.0;
	memset(f, -1, sizeof(f));
	dp(0, k + 1);
	printf("%.9lf\n", f[k + 1][0]);
	
	return 0;
}

递归+DP 103366 F Four Column Hanoi Tower

四层汉诺塔问题，实际上有有一个诡异的$O(1)$公式(我认为不需要掌握)，这里展示递推方法

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

int n;
int base = 10;
const int N = 1e4 + 10;
vector<int> p2[N], f[N];
// p2[i]: (2^i)-1

vector<int> Min(vector<int> A, vector<int> B){
	if (A.size() != B.size()) return (A.size() < B.size())?(A):(B);
	for (int i=A.size()-1;i>=0;--i){
		if (A[i] != B[i]) return (A[i] < B[i])?(A):(B);
	}
	return A;
}

vector<int> add(vector<int> A, vector<int> B){
    vector<int> res;
    if (A.size()<B.size()){
        return add(B,A);
    }
    int t = 0;
    for (int i=0;i<A.size();++i){
        t += A[i];
        if (i<B.size()){
            t += B[i];
        }
        res.push_back(t%base);
        t /= base;
    }
    if (t){
        res.push_back(1);
    }
    return res;
}

vector<int> mul(vector<int> &A, int b){
	vector<int> res;
	int t = 0;
	for (int i=0;i<A.size() | t;++i){
		if (i<A.size()) t += A[i]*b;
		res.push_back(t % base);
		t /= base;
	}
	while (res.size()>1 && res.back()==0){
		res.pop_back();
	}
	return res;
}

void init(){
	p2[0].push_back(0);
	vector<int> one;
	one.push_back(1);
	for (int i=1;i<200;++i){
		p2[i] = add(mul(p2[i-1], 2), one);
	}
	f[1].push_back(1);
	f[2].push_back(3);
	for (int i=3;i<=10000;++i){
		f[i] = p2[180];
		for (int j=1;j<min(150, i);++j){
			f[i] = Min(f[i], add(mul(f[i-j], 2), p2[j]) );
		}
	}
}

void solve(){
	scanf("%d", &n);
	
	for (int i=f[n].size()-1;i>=0;--i){
	    printf(i==0?"%d\n":"%d", f[n][i]);
	}
}

int main(){
	int t;
	scanf("%d", &t);
	init();
	while (t--){
		solve();
	}
	return 0;
}