与程序竞赛有关的数学知识点

内容较多,建议使用目录

容斥原理

集合的并集

设$A_1,A_2…,A_n$是有限集合,则
$|A_1\bigcup A_2\bigcup … \bigcup A_n|=\sum_{i=1}^n |A_i|-\sum_{i=1}^n \sum_{j>i}|A_i\bigcap A_j|-\sum_{i=1}^n \sum_{j>i} \sum_{k>j}|A_i\bigcap A_j\bigcap A_k|…+(-1)^n|A_1\bigcap A_2\bigcap…\bigcap A_n|$

Sylvester公式

给定集合N和具有性质i的集合$A_1,A_2…,A_n$,则
$|\overline A_1\bigcap \overline A_2\bigcap … \bigcap \overline A_n|=|N|-(\sum_{i=1}^n |A_i|-\sum_{i=1}^n \sum_{j>i}|A_i\bigcap A_j|…+(-1)^n|A_1\bigcap A_2\bigcap…\bigcap A_n|)$

[列题] 区间(S,E]中与n互质的元素个数

设$A_i$为n的第i个质因子$p_i$的倍数的集合且$A_i⊆(S,E],i=1,2,3…k$,那么$|A_i|=[\frac {E}{p_i}]-[\frac {S}{p_i}],|A_i\bigcap A_j|=[\frac {E}{p_ip_j}]-[\frac {S}{p_ip_j}]…$,然后套用上面公式即可解决问题了。关键是将这个公式表达出来,有没有发现一共有$2^n$项多项式,而且组合数量为奇数时负,组合数量为偶数时符号为正。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
LL POIAE(int n,int S,int E){
//质因子分解
int p[32],u=0;
for(int i=2;i*i<=n;i++)
if(n%i==0){
p[u++]=i;
while(n%i==0) n/=i;
}
if(n>1) p[u++]=n;

LL qua=0;
for(int i=1;i<1<<u;i++){
int cnt=0,com=1;
for(int k=0;k<u;k++)
if(i>>k&1) cnt++,com*=p[k];
qua+=cnt&1?E/com-S/com:-E/com+S/com;
}
return E-S-qua;
}


欧拉函数φ(n)

欧拉函数φ(n) 是小于等于n且与n互素的正整数的个数,假设$n=p_1^{α_1}p_2^{α_2}p_3^{α_3}…p_k^{α_k}$,则有
$$φ(n)=n(1-\frac{1}{p_1})(1-\frac{1}{p_2})(1-\frac{1}{p_3})…(1-\frac{1}{p_k})$$

证明:

设$A_i$为1到n之间$p_i$的倍数的集合,i=1,2,3,…,k,则
φ(n)=$|\overline A_1\bigcap \overline A_2\bigcap … \bigcap \overline A_k|$
=$|N|-(\sum_{i=1}^k |A_i|-\sum_{i=1}^k \sum_{j>i}|A_i\bigcap A_j|…+(-1)^k|A_1\bigcap A_2\bigcap…\bigcap A_k|)$
=$n-(\sum_{i=1}^k \frac{n}{p_i}-\sum_{i=1}^k \sum_{j>i}\frac{n}{p_i p_j}+\sum_{i=1}^k \sum_{j>i} \sum_{h>j}\frac{n}{p_i p_j p_h}…+(-1)^k\frac{n}{p_1 p_2 p_3…p_k})$
=$n(1-\frac{1}{p_1})(1-\frac{1}{p_2})(1-\frac{1}{p_3})…(1-\frac{1}{p_k})$

性质:

  • 欧拉函数是积性函数,即是说若m,n互质$φ(mn)=φ(m)φ(n)$ 。
  • 若n是质数p的k次幂,$φ(n)=φ(p^k)=p^k-p^{k-1}=(p-1)p^{k-1}$,因为除了p的倍数外,其他数都跟n互质。

欧拉函数φ(n) c++实现代码:

我们为了方便根据以上性质将其变形得到

$$ φ(n)=(p_1-1)p_1^{α_1-1}(p_1-1)p_2^{α_2-1}…(p_k-1)p_k^{α_k-1}$$

如此便可写出φ(n)的$o(\sqrt n )$ 的代码了。

1
2
3
4
5
6
7
8
9
10
11
LL Eular(int n){
LL ans=1;
for(int i=2;i*i<=n;i++)
if(n%i==0){
ans*=i-1;
n/=i;
while(n%i==0) n/=i,ans*=i;
}
if(n>1) ans*=n-1;
return ans;
}

欧拉函数值表筛法

1
2
3
4
5
6
7
int e[MAX_N];
void Euler(){
for(int i=0;i<MAX_N;i++) e[i]=i;
for(int i=2;i<MAX_N;i++)
if(e[i]==i)
for(int k=i;k<MAX_N;k+=i) e[k]=e[k]/i*(i-1);
}

[列1]区间[1,n]与n的最大公约数的和

$$\sum_{i=1}^n gcd(i,n)=\sum_{d|n}d\cdot φ(n/d) $$

[列2]

UVA11426 求 $\sum_{i=1}^n \sum_{j=i+1}^n gcd(i,j)$

1
2
3
4
5
6
7
8
9
10
11
12
13
int phi[MAX_N];LL res[MAX_N];
void solve(){//筛法的魅力
for(int i=1;i<MAX_N;i++) phi[i]=i;
for(int i=2;i<MAX_N;i++)
if(phi[i]==i)
for(int j=i;j<MAX_N;j+=i)
phi[j]=phi[j]/i*(i-1);
for(int i=1;i<MAX_N;i++)
for(int k=i;k<MAX_N;k+=i)
res[k]+=1LL*i*phi[k/i];
for(int i=1;i<MAX_N;i++)
res[i]+=res[i-1]-i;
}

[列3]

[1,n]与n互素的和
$$ \sum_{gcd(i,n)=1}^{i<n}i=n\cdot φ(n)/2 $$


莫比乌斯函数

莫比乌斯函数$μ(n)$,在狄利克雷卷积的乘法中与恒等函数互为逆元,$μ(1)=1$,对于无平方因子数$n=\prod_{i=1}^{t}p_i$有$\mu(n)=(-1)^t$,对于有平方因子数$n$有$μ(n)=0$。

$o(\sqrt n)$ 的板子

1
2
3
4
5
6
7
8
9
10
11
int mu(int n) {
int cnt = 0;
for(int i = 2; i * i <= n; ++i) {
if(n%i == 0) {
n /= i;++cnt;
if(n%i == 0) return 0;
}
}
if(n > 1) ++cnt;
return cnt&1?-1:1;
}


同余

两个整数 a,b,若它们除以正整数 m所得的余数相等,则称 a,b对于模 m同余,记作$ a\equiv b{\pmod {m}}$

其性质这里不赘述了,可以参看 维基


费马小定理

费马小定理是数论中的一个定理:假如a是一个整数,p是一个质数,那么 $a^{p}-a$ 是p的倍数,可以表示为
$$ a^{p}\equiv a{\pmod p}$$
如果a不是p的倍数,这个定理也可以写成
$$a^{p-1}\equiv 1{\pmod p}$$
[列题]计算 $2^{100}$ 除以13的余数

$ 2^{100}\equiv 2^{12\times 8+4}{\pmod {13}} \equiv (2^{12})^{8}\cdot 2^{4}{\pmod {13}}\equiv 1^{8}\cdot 16{\pmod {13}} \equiv 16{\pmod {13}}\equiv 3{\pmod {13}}$

则答案为3


欧拉定理

若$gcd(a,n)=1$,则$$a^{φ(n)} ≡ 1(mod \ n)$$


扩展欧几里得算法

给予二整数a、b,必存在有整数x、y使得ax + by = gcd(a,b)
对于这个问题,扩展欧几里得算法能够快速解出使上式成立(x,y),复杂度o(log(n))

1
2
3
4
5
6
7
8
9
int exGcd(int a,int b,int &x,int &y){
if (b==0){
x=1,y=0;
return a;
}
int q=exGcd(b,a%b,y,x);
y-=a/b*x;
return q;
}

通过上面代码可以得到一组特解,记为$(x_0,y_0)$,那么通解可以表示为

$$\begin{cases}
x = x_0 - {b}\cdot k \\
y = y_0 +{a}\cdot k \\
k∈\Bbb Z \\
\end{cases}$$

例题 POJ1061


模逆元

若$a\cdot b ≡ 1(mod \ m)$,称a和b对于模数m来说互为逆元

定理:a存在模m的乘法逆元的充要条件是gcd(a,m) = 1

证明:

  1. 充分性:由欧拉定理知$a^{φ(m)} ≡ 1(mod \ m)$,那么必有$b=a^{φ(m)-1}$ 使$a\cdot b ≡ 1(mod \ m)$
  2. 必要性:既然$a\cdot b ≡ 1(mod \ m)$,即有$a\cdot b-m\cdot k=1=k’\cdot gcd(a,m),k∈\Bbb Z,k’∈\Bbb Z$成立,即有$gcd(a,m)=1$

[例] 已知 a, m, 求b

将定义式改写成

$a\cdot b-1=m\cdot k,k∈\Bbb Z \\ ⇒ a\cdot b-m\cdot k=1,k∈\Bbb Z \\ $
利用扩展欧几里得算法不难求得一组解$(b’,k’)$ 使得

$a\cdot b’+m\cdot k’=gcd(a,m) $

结合定理得
若$gcd(a,m)≠ 1$,那么b无解
若$gcd(a,m)= 1$,那么$b=\frac{b’}{gcd(a,m)}=b’$

代码

那么也能写出代码了,若无答案返回-1,若有返回最小的正解

1
2
3
4
5
int mod_inverse(int a,int m){
int x,y;
if(exGcd(a,m,x,y)==1) return (m+x%m)%m;
return -1
}

o(n)得出[1,n]的(mod p)逆元

该算法要求所mod p为素数且n< p,记 $a$ 的逆元表示为 $a^{-1} $
设$k=p\%i$,假设以求得$k^{-1} \pmod p$,则

$i\cdot i^{-1}≡ k\cdot k^{-1}≡1\pmod p $
$⇒i\cdot i^{-1} ≡ (p\%i) \cdot (p\%i)^{-1} \pmod p $
$⇒i\cdot i^{-1} ≡ (p - [p/i]\cdot i) \cdot (p\%i)^{-1} \pmod p $
$⇒i\cdot i^{-1} ≡ ((i-1)\cdot p +p - [p/i]\cdot i) \cdot (p\%i)^{-1} \pmod p $
$⇒i\cdot i^{-1} ≡ i\cdot (p - [p/i]) \cdot (p\%i)^{-1} \pmod p $
由于p为素数,则$gcd(i,p)=1$,则
$i^{-1} ≡ (p - [p/i])\cdot (p\%i)^{-1} \pmod p$

取[0,p)中的解
$i^{-1} = (p - [p/i])\cdot (p\%i)^{-1} \% p$

代码为

1
2
3
4
int inv[MAXN];
inv[1] = 1;
for (int i = 2; i<MAXN; i++)
inv[i] = inv[MOD%i]*(MOD-MOD/i)%MOD;

孙子定理(中国剩余定理)

用现代数学的语言来说明的话,中国剩余定理给出了以下的一元线性同余方程组:
$$(S) :
\begin{cases}
x \equiv a_1 \pmod {m_1} \\
x \equiv a_2 \pmod {m_2} \\
\vdots \\
x \equiv a_n \pmod {m_n}
\end{cases}
$$
有解的判定条件,并用构造法给出了在有解情况下解的具体形式。

中国剩余定理说明:假设整数$m_1,m_2,…,m_n$其中任两数互质,则对任意的整数:$a_1,…,a_n$,方程组$(S)$有解,并且通解可以用如下方式构造得到:

  1. 设 ${ M=m_1 \times m_2\times \cdots \times m_n = \prod _{i=1}^nm_i}$ 是整数$m_1,m_2,…,m_n$的乘积,并设$ { M_{i}=M/m_{i},\;\;\forall i\in {1,2,\cdots ,n}} $,即$ { M_{i}} $是除了$m_i$以外的n − 1个整数的乘积。
  2. 设${ t_{i}=M_{i}^{-1}}$模$ m_{i}$的逆元: $t_i M_i \equiv 1 \pmod {m_i}, \; \; \forall i \in {1, 2, \cdots , n}$。
  3. 方程组$ (S)$的通解形式为:$x = a_1 t_1 M_1 + a_2 t_2 M_2 + \cdots + a_n t_n M_n + k M= k M + \sum_{i=1}^n a_i t_i M_i, \quad k \in \mathbb{Z}. $在模 $M$的意义下,方程组 $(S)$只有一个解:
    $$ x = \sum_{i=1}^n a_i t_i M_i.$$

证明

对于 ${ \forall j\in {1,2,\cdots ,n},\;j\neq i,\;\;\operatorname {gcd} (m_{i},m_{j})=1}$

考察乘积$ a_i t_i M_i $可知:

$ a_i t_i M_i \equiv a_i \cdot 1 \equiv a_i \pmod {m_i}, $
$\forall j \in {1, 2, \cdots , n}, \; j\neq i, \; \; a_i t_i M_i \equiv 0 \pmod {m_j}. $

所以 $x = a_1 t_1 M_1 + a_2 t_2 M_2 + \cdots + a_n t_n M_n$ 满足:

$\forall i \in {1, 2, \cdots , n}, \; \; x = a_i t_i M_i + \sum_{j \neq i} a_j t_j M_j \equiv a_i + \sum_{j \neq i} 0 \equiv a_i \pmod {m_i}. $

到此以得证$(S)$的一个解$ x = \sum_{i=1}^n a_i t_i M_i.$

通解

假设 $x_1$和$x_2$都是方程组 $(S)$的解,那么:

$ \forall i \in {1, 2, \cdots , n}, \; \; x_1 - x_2 \equiv 0 \pmod {m_i} .$,即 $\forall m_i|x_1-x_2$

即 $\prod_{i=1}^n m_i|x_1-x_2$,即 $M|x_1-x_2$

到此说明方程组 $(S)$ 的任何两个解之间必然相差$M$的整数倍,则通解可表示为
$${k M + \sum_{i=1}^n a_i t_i M_i \; ; \quad k \in \mathbb{Z} }. $$

代码

1
2
3
4
5
6
7
8
9
10
LL CRT(int *m,int *a,int n){
LL M=1,ans=0,x,t_i;
for (int i=1;i<=n;i++) M*=m[i];
for (int i=1;i<=n;i++){
LL M_i=M/m[i];
exGcd(m[i],M_i,x,t_i); //求 M_i 模 m[i] 的逆元 t_i
ans=(ans+a[i]*t_i*M_i)%M;
}
return ans;
}

矩阵

数学上,一个m×n的矩阵是一个由m行n列元素排列成的矩形阵列。矩阵里的元素可以是数字、符号或数学式。

大小相同(行数列数都相同)的矩阵之间可以相互加减,具体是对每个位置上的元素做加减法。矩阵的乘法则较为复杂。两个矩阵可以相乘,当且仅当第一个矩阵的列数等于第二个矩阵的行数。矩阵的乘法满足结合律和分配律,但不满足交换律。

矩阵乘法

矩阵相乘最重要的方法是一般矩阵乘积。它只有在第一个矩阵的列数(column)和第二个矩阵的行数(row)相同时才有定义。一般单指矩阵乘积时,指的便是一般矩阵乘积。若A为 ${ m\times n} $矩阵,B为$ { n\times p} $矩阵,则他们的乘积AB(有时记做A · B)会是一个 ${ m\times p} $ 矩阵。其乘积矩阵的元素如下面式子得出:
$$ AB_{ij}=S_{r=1}^n a_{ir} b_{rj} = a_{i1} b_{1j} + a_{i2} b_{2j} + \cdots + a_{in}b_{nj} $$

矩阵快速幂

类比快速幂不难的到矩阵快速幂,若有矩阵a,则求 $a^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
int N;
struct matrix{LL m[MAX_N][MAX_N];};
matrix multi(matrix a,matrix b){
matrix tmp;
for(int i=1;i<=N;i++)
for(int j=1;j<=N;j++){
tmp.m[i][j]=0;
for(int k=1;k<=N;k++)
tmp.m[i][j]+=a.m[i][k]*b.m[k][j];
tmp.m[i][j]%=MOD;
}
return tmp;
}
matrix fast_mod(matrix a,int n){
matrix res;
for(int i=1;i<=N;i++)
for(int j=1;j<=N;j++)
res.m[i][j]=(i==j);
while(n){
if(n&1) res=multi(res,a);
a=multi(a,a);
n>>=1;
}
return res;
}

线性递推关系与矩阵乘法

设数列 ${h_n}$ 满足 $k$ 阶常系数线性递推关系:
$$h_n = C_1h_{n−1} + C_2h_{n−2} + C_3h_{n−3} +···+ C_kh_{n−k}+b_n$$
$b_n$ 可以为常数,也可以是关于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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#include<stdio.h>
#include<string.h>
#include<queue>
#include<vector>
#define INF 0x3f3f3f3f
#define MAX_N 10
#define MOD 1000000007
using namespace std;
typedef long long LL;
int N;LL b_n=0,C[MAX_N],h[MAX_N];
struct matrix{LL m[MAX_N][MAX_N];};
matrix multi(matrix a,matrix b){
matrix tmp;
for(int i=1;i<=N;i++)
for(int j=1;j<=N;j++){
tmp.m[i][j]=0;
for(int k=1;k<=N;k++)
tmp.m[i][j]+=a.m[i][k]*b.m[k][j];
tmp.m[i][j]%=MOD;
}
return tmp;
}
matrix fast_mod(matrix a,int n){
matrix res;
for(int i=1;i<=N;i++)
for(int j=1;j<=N;j++)
res.m[i][j]=(i==j);
while(n){
if(n&1) res=multi(res,a);
a=multi(a,a);
n>>=1;
}
return res;
}
void init(matrix &res,matrix &H){
for(int i=1;i<=N;i++) res.m[1][i]=C[i];
res.m[1][N]=b_n;
for(int i=2;i<=N;i++)
for(int j=1;j<=N;j++)
res.m[i][j]=(i==j+1);
res.m[N][N-1]=0;res.m[N][N]=1;

for(int i=1;i<=N;i++){
H.m[i][1]=h[N-i];
for(int j=2;j<=N;j++) H.m[i][j]=0;
}
H.m[N][1]=1;
}
void slove(int k,int n){
matrix res,H;
init(res,H);

res=fast_mod(res,n-k+1);
res=multi(res,H);
LL ans=res.m[1][1];

printf("%lld\n",ans);
}

int main()
{
int k,n;
while(scanf("%d%d",&n,&k)!=EOF){
for(int i=1;i<=k;i++) scanf("%lld",&C[i]);
for(int i=1;i<=k;i++) scanf("%lld",&h[i]);
N=k+1;
slove(k,n);
}
return 0;
}

例题 51NOD 1126

###《挑战》上的矩阵快速幂模板
以斐波那契数列第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
#include<stdio.h>
#include<vector>
using namespace std;

typedef long long LL;
typedef vector<LL> vec;
typedef vector<vec> mat;
const LL MOD = 1000000009;

mat mul(mat &A,mat &B){
mat C(A.size(),vec(B[0].size()));
for(int i=0;i<A.size();++i){
for(int k=0;k<B.size();++k){
for(int j=0;j<B[0].size();++j)
C[i][j]=(C[i][j]+A[i][k]*B[k][j])%MOD;
}
}
return C;
}

mat pow(mat A,LL n){
mat B(A.size(),vec(A.size()));
for(int i=0;i<A.size();++i) B[i][i]=1;
while(n){
if(n&1) B=mul(B,A);
A=mul(A,A);
n>>=1;
}
return B;
}

int main()
{
mat A(2,vec(2));
A[0][0]=A[0][1]=A[1][0]=1;
A[1][1]=0;
LL n;
scanf("%lld",&n);
A=pow(A,n);
printf("%lld\n",A[1][0]);
return 0;
}


##其他

###约数个数
仍何一个整数n都可以表示成$n=p_1^{a_1}p_2^{a_2}…p_k^{a_k}$,那么,n的约数个数
$$\sum_{d|n} 1=(a_1+1)(a_2+1)…(a_k+1)$$
例题 LOJ1341

###a的k次幂前n位数

用该方法要注意精度

$a^k=10^{\lg{a^k}}=10^{k\cdot \lg a}=10^{[k\cdot \lg a]+{k\cdot \lg a}}$

前n位:[10^{n-1}\cdot 10^{\{k\cdot \lg a\}}

例题 LOJ1282

###定积分
Simpson`s 3/8 rule

1
2
3
4
5
6
7
8
9
10
11
12
13
14
const double eps=1e-8;
inline double f(double x){
return x*x;//被积函数
}
inline double getAppr(double a,double b){//Simpson`s 3/8 rule
return (b-a)*(f(a)+3*f((2*a+b)/3)+3*f((a+2*b)/3)+f(b))/8.0;
}
double simpson(double l,double r){//积分区间[l,r]
double sum=getAppr(l,r);
double mid=(l+r)/2;
double suml=getAppr(l,mid);
double sumr=getAppr(mid,r);
return fabs(sum-suml-sumr)<eps?sum:simpson(l,mid)+simpson(mid,r);
}

例题 地狱飞龙


##线性筛法筛素数、莫比乌斯函数、欧拉函数

线性筛法(欧拉筛法)可以在 O(n) 的时间内获得 [1,n] 的所有素数。算法保证每个合数都会被它的最小质因子筛掉,所以复杂度是线性的。同时,我们可以利用这一特性,结合积性函数的性质,在 O(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
bool isNotPrime[MAXN + 1];
int mu[MAXN + 1], phi[MAXN + 1], primes[MAXN + 1], cnt;
inline void euler()
{
isNotPrime[0] = isNotPrime[1] = true;
mu[1] = 1;
phi[1] = 1;
for (int i = 2; i <= MAXN; i++) {
if (!isNotPrime[i]) {
primes[++cnt] = i;
mu[i] = -1;
phi[i] = 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;
phi[t] = phi[i] * primes[j];
break;
} else {
mu[t] = -mu[i];
phi[t] = phi[i] * (primes[j] - 1);
}
}
}
}


本文章不定期更新

如有谬误,恳请指正