数学提高

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


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

常用数学技巧和公式

1. 常用公式

完全立方公式

(a+b)3=a3+3a2b+3ab2+b3=a3+b3+3ab(a+b)(a+b)^3=a^3+3a^2b+3ab^2+b^3=a^3+b^3+3ab(a+b)

a3+b3=(a+b)33ab(a+b)=(a+b)(a2+b2ab)a^3+b^3=(a+b)^3-3ab(a+b)=(a+b)(a^2+b^2-ab)

(ab)3=a33a2b+3ab2b3=a3b3+3ab(ba)(a-b)^3=a^3-3a^2b+3ab^2-b^3=a^3-b^3+3ab(b-a)

a3b3=(ab)3+3ab(ab)=(ab)(a2+b2+ab)a^3-b^3=(a-b)^3+3ab(a-b)=(a-b)(a^2+b^2+ab)

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

1. 拆数

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

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

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

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

(x1)3+(x1)3x3x3=6x(x-1)^3+(x-1)^3-x^3-x^3=6x

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

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

假设x=2n+1x=2n+1,显然2n+1=(n+1)2n22n+1=(n+1)^2-n^2

随机化

mt19937生成随机数

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

函数用法:

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

完整代码

1
2
3
4
5
6
7
8
9
10
11
12
#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(fn,fm)\gcd(f_n, f_m)

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#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;
}

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
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) 应用范围:递推式操作

一般表现为如下形式:某个行向量Xn=[an  an1 ..ank+1]X_n=[a_n~~ a_{n-1}~..a_{n-k+1}],边界X1=[ak1  ak2 ..a0]X_1=[a_{k-1}~~a_{k-2}~..a_0],可以找到某个矩阵AA,满足:Xn=Xn1×AX_n=X_{n-1}\times A,则可以得到Xn=X1×An1X_n=X_1\times A^{n-1}。注意:一定要保证AA中没有变量。

由于矩阵具有结合律,我们可以先求出B=An1modPB=A^{n-1}\mod P,然后再计算X1×BX_1\times B,即可求出XnX_nXnX_n的第一个元素就是ana_n

(3) 时间复杂度:O(logn)O(\log n) 但是大常数ORZ

(4) 例题+模板

① AcWing1303 斐波那契前n项和

现有Fibonacci数列如下:f1=1,f2=1,fn=fn2+fn1f_1=1,f_2=1,f_n=f_{n-2}+f_{n-1}

定义Sn=i=1nfiS_n=\sum\limits _{i=1}^n f_i,对于给定的n,mn,m,求SnmodmS_n \mod m

数据范围:1n2×109;1m109+101\leq n \leq 2\times 10^9;1\leq m \leq 10^9 + 10

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

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

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

那么[fn,fn+1,Sn]×[0  1  01  1  10  0  1]=[fn+1,fn+2,Sn+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=[0  1  01  1  10  0  1]A= \begin{bmatrix} 0~~1~~0 \\ 1~~1~~1 \\ 0~~0~~1 \end{bmatrix},求出An1modmA^{n-1}\mod m后,计算X1×An1modmX_1\times A^{n-1}\mod m即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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)=(f1+2f2+..+nfn)modmT(n)=(f_1+2f_2+..+nf_n)\mod m的值,数据范围:1n,m23111\leq n,m\leq 2^{31}-1

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

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

{nSnTn=(n1)f1+(n2)f2+..+fn1(n+1)Sn+1Tn+1=nf1+(n1)f2+..+2fn1+fn\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}

定义Pn=nSnTnP_n=nS_n-T_n可得Pn+1Pn=SnP_{n+1}-P_n=S_n,而Tn=nSnPnT_n=nS_n-P_n

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

[fn,fn+1,Sn,Pn]×[0  1  0  01  1  1  00  0  1  10  0  0  1]=[fn+1,fn+2,Sn+1,Pn+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=[0  1  0  01  1  1  00  0  1  10  0  0  1]A=\begin{bmatrix} 0~~1~~0~~0\\ 1~~1~~1~~0\\ 0~~0~~1~~1\\ 0~~0~~0~~1 \end{bmatrix},然后算出Pn,SnP_n,S_n,则Tn=nSnPnT_n=nS_n-P_n,本题结束。

③ P7453 大魔法师

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

④ Loj6208 树上询问

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

⑤ AcWing1305 GT考试

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

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

假如有090-9一共1010个数, 我们在第i+1i+1个位置填cc, 就可以利用kmpkmp求出来此时匹配到了第jj个位置
也就是说, 当前已经填好了第ii个数, 准备在第i+1i+1的位置填cc的时候, 可以转移到jj
f[i][k]=f[i1][0]×a[0][k]+f[i1][1]×a[1][k]+..+f[m1][k]×a[m1][k]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>ji->j的贡献, 就得到了我们的转移矩阵
Xi=Xi1×AX_{i}= X_{i-1} \times A

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
#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序列统计的形式简单讲下。

image-20220717181814104

约定:题目中的 N,L,RN,L,Rn,l,rn,l,r 表示。

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

1. 映射法

image-20220717182521922

image-20220717182534272

2. 隔板法

image-20220717182618345

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

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#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

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

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

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

同样的,假如我们已经处理好i+1i+1nn的答案,rir_i也能不断的向右扩展。

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

那么完整代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
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

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

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
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×mn\times m的网格,里面有最多A to ZA~to ~Z个孩子,king要把这个网格分成若干矩形,要满足:

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

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

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

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
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(n)O(\sqrt n)的复杂度快速计算形如i=1nf(i)g(ni)\sum \limits_{i=1}^n f(i)g(\lfloor \frac n i \rfloor)的含有除法向下取整的和式,它的核心思想是将ni\lfloor \frac n i \rfloor相同的数打包同时计算,主要利用了Fubini定理。

证明

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

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

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

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

2. 证明算法的正确性

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

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

t>l,t<r\because t>l,t<r

ntnl,ntnr\therefore \lfloor \frac n t \rfloor \geq \lfloor \frac n l \rfloor, \lfloor \frac n t \rfloor \leq \lfloor \frac n r \rfloor ,即nt=nl\lfloor \frac n t \rfloor=\lfloor \frac n l \rfloor,矛盾,故得证。

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

显然,C=nlC=\lfloor \frac n l \rfloor,下面就要证明r=nnl=nCr=\lfloor \frac{n}{\lfloor \frac n l \rfloor} \rfloor=\lfloor \frac n C \rfloor

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

  • pp满足ii的性质,即np=C\lfloor \frac n p \rfloor=C

p=nkCp=\frac{n-k}{C},则np=C×nnkC×nn\frac n p=\frac{C\times n}{n-k}\geq \frac {C\times n}{n},即npC\frac n p \geq C

nmodp=k,k<pn\mod p=k,k<p,则nk<p×(C+1)n-k<p\times (C+1),即np<C+1\lfloor \frac n p \rfloor<C+1

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

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

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

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

p×C+C=ng\therefore p\times C + C=n-g,又p×C=nkp\times C=n-k

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

综上,nC\lfloor \frac n C \rfloor即为令ni=C\lfloor \frac n i \rfloor=C的最大的数

模板

1
2
3
4
5
6
7
8
9
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,kn,k,计算G(n,k)=i=1nkmodiG(n,k)=\sum\limits _{i=1}^n k \bmod i

数据范围:1n,k1091\leq n,k\leq 10^9

(懒得自己写题解了,直接cpy洛谷题解orz)由题意得:ans=i=1nkmodians=∑_{i=1}^nk\bmod i

我们知道,amodb=ab×aba \bmod b = a - b \times \lfloor \frac{a}{b} \rfloor

因此,ans=i=1nki×ki=nki=1ni×kians = \sum_{i=1}^{n} k - i \times \lfloor \frac{k}{i} \rfloor = nk - \sum_{i=1}^{n} i \times \lfloor \frac{k}{i} \rfloor

首先枚举块的左边界 ll,并根据左边界和$ k$ 计算出右边界 rr

t=klt = \lfloor \frac{k}{l} \rfloor,分两种情况讨论:

  • t0t \neq 0,则 r=min(kt,n)r = \min (\lfloor \frac{k}{t} \rfloor , n)
  • t=0t = 0,则 r=nr = n

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

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
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)f(n),称为一个数论函数,又叫算术函数。:

**举例:potp(n),μ(n),φ(n),π(n)pot _p (n),\mu(n),\varphi(n),\pi(n)**都是数论函数,n(i)=nn(i)=n:一直返回同一个常数nn

加性函数与积性函数

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

则称ff加性的,这时定义给出f(1)=0f(1)=0

  • m,n\forall m,n,满足gcd(m,n)=1gcd(m,n)=1,有f(m×n)=f(m)×f(n)f(m\times n)=f(m)\times f(n)则称ff积性的,这时定义给出f(1)=1f(1)=1

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

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

nn的唯一分解为n=i=1kpiai=p1a1p2a2..pkakn= \prod_{i=1}^kp_i^{a_i} = p_1^{a_1}p_2^{a_2}..p_k^{a_k},那么

  • f(1)=1f(1)=1f(n)=i=1kf(piai)f(n)f(n)=\prod_{i=1}^k f(p_i^{a^i}) \Leftrightarrow f(n)积性函数(充要条件)

  • f(1)=1f(1)=1fa1(p1)×fa2(p2)×..×fak(pk)f(n)f^{a_1}(p_1)\times f^{a_2}(p_2)\times..\times f^{a_k}(p_k) \Leftrightarrow f(n)完全积性函数(充要条件)

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

举例:常见的积性函数

单位元函数e(n)=[n=1]e(n)=[n=1],即e(n)={1n=10n1e(n)=\begin{cases} 1 & n=1\\ 0 & n\neq 1 \end{cases}

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

幂函数idk(n)=nkid^k(n)=n^k,完全积性

单位函数id(n)=nid(n)=n,完全积性,相当于id1id^1

恒等函数I(n)=1I(n)=1(也就是常数值函数,只是刚好这个常数等于1),相当于id0id^0

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

除数函数σk(n)=dndk\sigma_k(n)=\sum\limits _{d|n}d^k,表示nn的约数的kk次幂和

约数和函数σ(n)=σ1(n)=dnd\sigma(n)=\sigma_1(n)=\sum\limits _{d|n}d,表示nn的约数之和

约数个数函数τ(n)=σ0(n)=dn1\tau(n)=\sigma_0(n)=\sum\limits _{d|n}1,表示nn的约数的个数,一般也写作d(n)d(n)

莫比乌斯函数μ(n)\mu(n)

生成函数

image-20220717185858084

狄利克雷卷积

定义

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

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

h(n)=dnf(d)g(nd)h(n)=\sum\limits _{d|n}f(d)g(\frac n d),或者说h(n)=xy=nf(x)g(y)h(n)=\sum\limits _{xy=n}f(x)g(y)

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

2. 狄利克雷逆(Dirichlet inverse)

fg=ef*g=e,则称ggff的狄利克雷逆,记作f1f^{-1}

计算狄利克雷逆

首先,当n=1n=1

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

(ff1)(1)=e(1)=1(f*f^{-1})(1)=e(1)=1,所以f(1)f1(1)=1f(1)f^{-1}(1)=1

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

这说明了,f1f^{-1}存在的必要条件f(1)0f(1)\neq 0

n>1n>1时,有

(ff1)(n)=dnf(d)f1(nd)=f(1)f1(n)+dn,d>1f(d)f1(nd)(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)

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

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

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

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

f1(n)={1f(1)n=11f(1)dn,d>1f(d)f1(nd)n>1f^{-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}

性质

交换律

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

结合律

image-20220717183440979

分配律

image-20220717183452655

f,gf,g为积性函数,则h=fgh=f*g也是积性函数

image-20220717183553371

g,h=fgg, h=f*g是积性函数,则ff也是积性函数

image-20220717183628289

数论函数的卷积关系

0. 对于任意数论函数ff和恒等函数II,有

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

1. 幂函数与除数函数

(idkI)(n)=dnidk(n)=nk=σk(n)(id^k*I)(n)=\sum\limits _{d|n}id^k(n)=\sum\limits n^k=\sigma_k(n)

idkI=σkid^k * I=\sigma_k

2. 欧拉函数与恒等函数

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

3. 恒等函数与恒等函数

II=dI *I=d(在莫比乌斯函数里证明过了)

莫比乌斯反演

莫比乌斯函数

1. 定义

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

μ(n)={1n=10i,ai2(1)ki,ai=1\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

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

证明:设ddnn的约数,则d=i=1kpibid=\prod_{i=1}^kp_i^{b_i},其中0biai0\leq b_i\leq a_i

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

dnμ(d)=i=0kCki×(1)i\sum\limits _{d|n}\mu(d)=\sum\limits _{i=0}^kC_k^i\times(-1)^i,又(ab)k=i=0kCkiakbki(a-b)^k=\sum\limits_{i=0}^k C_k^i a^kb^{k-i}

(11)k=i=0kCki×(1)k(1-1)^k=\sum\limits _{i=0}^kC_k^i\times (-1)^k,故dnμ(d)=0k=0\sum\limits _{d|n}\mu(d)=0^k=0

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

(1) μI=e\mu * I = e

证明:设n=i=1kpiai,n=i=1kpin=\prod_{i=1}^kp_i^{a_i}, n'=\prod_{i=1}^kp_i

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

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

因此,μ\muII的狄利克雷逆。

(2) μid=φ\mu * id = \varphi

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

(3) μd=I\mu * d=I

证明:(II)(n)=dnI(d)=dn1=d(n)(I*I)(n)=\sum\limits _{d|n}I(d)=\sum\limits _{d|n}1=d(n)

d=II\therefore d=I*I,又μ=I1\mu=I^{-1}

μd=I\therefore \mu * d=I

4. 线性筛法求莫比乌斯函数
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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(n)=fI=dnf(d)f(n)=dnμ(d)×F(nd)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)=fIF(n)=f*I,则f=FI1=Fμf=F*I^{-1}=F*\mu

f(n)=dnμ(d)×F(nd)f(n)=\sum\limits_{d|n}\mu(d)\times F(\frac n d)

证明(利用性质1+二重积分交换次序的思想):

dnμ(d)×F(nd)=dnμ(d)×indf(i)=inf(i)dniμ(d)\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)

ii能取到所有dd可以取到的取值,这样反过来看,把ii提到前面)

又当且仅当n=in=i时,dniμ(d)=1\sum\limits_{d|\frac{n}i}\mu(d)=1,因此dnμ(d)×F(nd)=f(n)\sum\limits_{d|n}\mu(d)\times F(\frac n d)=f(n)

倍数形式:

  • F(n)=nNf(N)f(n)=nNF(N)μ(Nn)F(n)=\sum\limits_ {n|N}f(N) \Leftrightarrow f(n)=\sum\limits_{n|N}F(N)\mu(\frac N n),(枚举NNnn的所有倍数,N[n,+)N\in[n,+\infin)

证明:nNF(N)μ(Nn)=nNμ(Nn)Nif(i)\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=Nnd=\frac N n,则N=dnN=dn,则dnidn|i,即dind|\frac i n

因此nNμ(Nn)Nif(i)=dinμ(d)Nif(i)\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=in=i时,dniμ(d)=1\sum\limits_{d|\frac{n}i}\mu(d)=1,因此f(n)=nNF(N)μ(Nn)f(n)=\sum\limits_{n|N}F(N)\mu(\frac N n)

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

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

p2522 Problem b

image-20220717183649103

数据范围:1n,k5×104;1ab5×104;1cd5×1041\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,c1]A[1,a1][1,d]+A[1,a1][1,c1]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]

image-20220717183727183

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

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

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

f(k)=kNμ(Nk)ak×bkf(k)=\sum\limits _{k|N}\mu(\frac N k)\lfloor \frac a k \rfloor\times \lfloor \frac b k \rfloor,设 t=Nkt=\frac{N}{k},显然枚举 tt的结果为1,2,..,1,2,..,这样的整数,N=tkN=tk

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
#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)d(x)xx的约数个数,给定TTn,mn,m,求i=1Nj=1Md(i×j)\sum\limits _{i=1}^N \sum\limits_{j=1}^M d(i\times j)

数据范围:1N,M,T5×1041\leq N,M,T\leq 5\times 10^4

i=1Nj=1Md(i×j)=i=1Nj=1Mxiyj[(x,y)=1]\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=i=1kpiai,j=i=1kpibii=\prod_{i=1}^k p_i^{a_i},j=\prod_{i=1}^k p_i^{b_i}0ai,bi0\leq a_i,b_i

i×j=i=1kpiai+bii\times j=\prod_{i=1}^k p_i^{a_i+b_i}d(i×j)=i=1k(ai+bi+1)d(i\times j)=\prod_{i=1}^k(a_i+b_i+1)

即从ii中选出约数xxjj中选出约数yy,对于p1p_1而言,若要求(x,y)=1(x,y)=1

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

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

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

f(n)=i=1Nj=1Mxiyj[(x,y)=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(1)就是答案。

F(n)=i=1Nj=1Mxiyj[n(x,y)]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)=ndf(d)F(n)=\sum\limits _{n|d}f(d)

f(n)=ndμ(dn)F(d)f(n)=\sum\limits _{n|d}\mu(\frac d n)F(d)T=min(N,M)T=\min(N,M),则f(1)=d=1Tμ(d)F(d)f(1)=\sum\limits _{d=1}^T\mu(d)F(d)

  • 再化简FF

F(n)=i=1Nj=1Mxiyj[n(x,y)]=x=1Ny=1MNxMy[n(x,y)]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)]

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

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

x=xn,y=ynx'=\lfloor \frac x n \rfloor,y'=\lfloor \frac y n \rfloor

F(n)=x=1Ny=1MNxMy[n(x,y)]=x=1Nny=1MnNnxMnyF(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=Nn,M=MnN'=\lfloor \frac N n \rfloor,M'=\lfloor \frac M n \rfloor

F(n)=x=1Ny=1MNxMy=(x=1NNx)×(y=1MMy)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)=i=1nni)h(n)=\sum\limits_{i=1}^{n} \lfloor \frac {n} {i} \rfloor),也就是标准整数分块,则F(n)=h(N)×h(M)F(n)=h(N')\times h(M')

  • 再求f(1)f(1)

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
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={p1,p2,..,pk}S=\{p_1,p_2,..,p_k\},只由这些质数相乘得到的数我们成为丑数(Humble/Ugly Numbers)。习惯上,我们认为第一个丑数是1,但是也可能不是,所以看清题意。记第nn大 的丑数为h[n]h[n]

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

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#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×k)O(n\times k)注意,这里的下标是从0开始的,并且,INF要足够足够大!!!!

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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. 正常的快速幂
1
2
3
4
5
6
7
8
9
10
// 求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公式:AbAbmodφ(m)(modm)A^b\equiv A^{b\mod \varphi(m)} (\mod m),即可。因为经常用到高精度,所以在这里放一个Python的代码吧

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
# -*- 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. 十进制快速幂

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

比如,3405=(31)5×(310)0+(3100)43^{405}=(3^1)^5\times(3^{10})^0+(3^{100})^4

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#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(logb)O(\log b)
1
2
3
4
5
6
7
8
9
10
// 求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)O(1)

应用范围:计算a×bmodm;a,bm1018a\times b \bmod m;a,b\leq m\leq 10^{18}

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

image-20220717183820362

1
2
3
4
5
6
7
8
9
10
#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

给定NN,要求构造出x13+x23+...+xk3=Nx_1^3+x_2^3+...+x_k^3=Nk5k\leq 5

数据范围:1N10100000;xi101100001\leq N\leq 10^{100000};|x_i|\leq 10^{110000}

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

认真的说,这个更像完全平方公式的拓展:(x+1)2(x1)2x2x2=4x(x+1)^2-(x-1)^2-x^2-x^2=4x