数论
一、质数
质数的定义：在大于1的整数中，只包含1和它本身两个约数的数称为质数，也称作素数。
1、质数的判定
（1）试除法：根据定义来判断
最朴素的方法，时间复杂度O(n)，代码如下：
bool is_prime(int n)
{
if(n < 2)
return false;
for(int i = 2;i < n;i++)
if(n % i == 0)
return false;
return true;
}
对此方法的优化：我们可以发现，一个数所有的约数都是成对出现的，所以我们可以假设 a × b = m a\times b=m a×b=m那么我们可以只找小于 x × x = m x \times x=m x×x=m的部分，每次都找一对即可，这样只需要找到根号n就行，时间复杂度为 O n O\sqrt n On，代码为：
bool is_prime(int n)
{
if(n < 2)
return false;
for(int i = 2;i <= n/i;i++) //这里不用sqrt是为了更快一点，每次调用sqrt函数会慢
if(n % i == 0)
return false;
return true;
}
（2）分解质因数法：将n分解为他的所有质因数
我们将n分解质因数，使得
n
=
β
1
α
1
×
β
2
α
2
×
.
.
.
×
β
k
α
k
n = \beta_1^{\alpha_1}\times \beta_2^{\alpha_2}\times...\times \beta_k^{\alpha_k}
n=β1α1×β2α2×...×βkαk
首先我们可以知道对于一个数n中最多只有一个大于根号n的质因子，（如果有两个那么相乘都已经大于n了），这样我们就可以只找小于根号n的质因子，最后再判断是否有大于根号n的，时间复杂度为
O
(
n
)
O(\sqrt n)
O(n)，分解输出质因子和幂次的代码如下：
bool is_prime(int n)
{
for(int i = 2;i <= n/i;i++)
{
if(n % i == 0) //i是质因数
{
int s = 0;
while(n % i == 0) //循环计算次数
{
n /= i;
s++;
}
cout << i << " " << s << endl; //i是质因数，s是幂次数
}
}
if(n > 1) //如果n还大于1说明有大于根号n的质因数
cout << n << " " << 1 << endl; //就是n
}
2、质数筛
（1）埃氏筛：将小于等于n的所有质数筛出来
我们循环每一个数，将他的所有倍数筛掉（一定不是质数），则剩下的一定都是质数，时间复杂度： O ( n log log n ) O(n \log\log n) O(nloglogn)可近似为 O ( n ) O(n) O(n)，代码如下：
bool vis[MAXN];
int primes[MAXN],cnt;
bool get_primes(int n)
{
for(int i = 2;i <= n;i++)
{
if(!vis[i])
{
primes[cnt++] = n;
for(int j = i+i;j <= n;j += i)
vis[j] = 1;
}
}
}
（2）欧拉筛（线性筛）
我们看埃氏筛，对于某个合数，我们可能将他筛了很多次，（循环到他的每一个质因数时都筛到他了），因此效率肯定低，那么就需要找一种办法来解决，就是欧拉筛。
欧拉筛又称线性筛，他是在线性（
O
(
n
)
O(n)
O(n)）的情况下筛完所有素数。其核心思想是只用每个合数的最小质因子去筛掉该合数，这样会把每个合数都筛掉（每个合数都有其最小质因子），保证了正确性，同时对每个数都只访问了一遍，是线性的。
int primes[MAXN],cnt;
bool vis[N];
void get_primes(int n)
{
for(int i = 2;i <= n;i++)
{
if(!vis[i]) //如果i没有被筛掉，那么他就是质数
primes[cnt++] = i;
for(int j = 0;primes[j] <= n/i;j++) //从小到大循环所有质数
{
vis[primes[j]*i] = 1; //筛掉primes[j]*i
if(i % primes[j] == 0) //线性的关键语句
break;
}
}
}
该代码的关键就是if(i % primes[j] == 0) break;
①对于i % primes[j] == 0的情况时，因为是从小到大枚举的质数，所以primes[j]一定是i*primes[j]的最小质因数
②对于i % primes[j] != 0的情况时，因为是从小到大枚举的质数primes[]，还没有碰到i % primes[j] == 0的情况，因此当前的primes[j]一定比i的最小质因子小，因此primes[j]一定是i*primes[j]的最小质因数
综上所述，任何情况时，primes[j]都保证了一定是i*primes[j]的最小质因数，而根据代码vis[primes[j]*i] = 1;我们每次筛掉的合数就是i*primes[j]，这样就保证了我们对于每个合数一定是由他的最小质因数筛掉的，就保证了线性的复杂度；而同时，我们又保证了循环每个素数作为最小质因数去筛，就会把所有的合数筛掉，也保证了筛法的正确性。
二、约数
首先我们要知道，对于任意一个数n，可以将他质因数分解为： n = β 1 α 1 × β 2 α 2 × . . . × β k α k n = \beta_1^{\alpha_1}\times \beta_2^{\alpha_2}\times...\times \beta_k^{\alpha_k} n=β1α1×β2α2×...×βkαk，其中 β i \beta _i βi是质因数。（也可以是质数，非质因数的质数的幂次一定是0）
1、约数个数
约数的个数有：
(
α
1
+
1
)
(
α
2
+
1
)
.
.
.
(
α
k
+
1
)
(\alpha_1+1)(\alpha_2+1)...(\alpha_k+1)
(α1+1)(α2+1)...(αk+1)个。
因为根据组合数学，对于每一个质因数
β
i
\beta _i
βi他的指数可以是
0
,
1...
α
i
0,1...\alpha_i
0,1...αi，共
α
i
+
1
\alpha_i+1
αi+1项，因此总共的排列组合结果有
(
α
1
+
1
)
(
α
2
+
1
)
.
.
.
(
α
k
+
1
)
(\alpha_1+1)(\alpha_2+1)...(\alpha_k+1)
(α1+1)(α2+1)...(αk+1)项，即约数个数。
2、约数之和
约数之和为：
(
β
1
0
+
β
1
1
+
.
.
.
+
β
1
α
1
)
(
β
2
0
+
β
2
1
+
.
.
.
+
β
2
α
2
)
.
.
.
(
β
k
0
+
β
k
1
+
.
.
.
+
β
k
α
k
)
(\beta_1^0+\beta_1^1+...+\beta_1^{\alpha_1})(\beta_2^0+\beta_2^1+...+\beta_2^{\alpha_2})...(\beta_k^0+\beta_k^1+...+\beta_k^{\alpha_k})
(β10+β11+...+β1α1)(β20+β21+...+β2α2)...(βk0+βk1+...+βkαk)。
我们可以将上式整个拆开会得到很多项相加的形式，根据排列组合定理，一共有
(
α
1
+
1
)
(
α
2
+
1
)
.
.
.
(
α
k
+
1
)
(\alpha_1+1)(\alpha_2+1)...(\alpha_k+1)
(α1+1)(α2+1)...(αk+1)项，而每一项正好就是n的一个约数，则其就是所有约数相加，即约数之和。
因此约数个数和约数之和其实就是分解质因数，分解质因数的方法在上文有提。
三、欧拉定理
1、欧拉函数
（1）定义
欧拉函数记作
φ
(
n
)
\varphi (n)
φ(n)，表示小于n的正整数中与n互质的个数。互质指最大公因数为1。
（2）公式
那这个函数怎么求呢，其公式为：若n能因式分解为
n
=
β
1
α
1
×
β
2
α
2
×
.
.
.
×
β
k
α
k
n = \beta_1^{\alpha_1}\times \beta_2^{\alpha_2}\times...\times \beta_k^{\alpha_k}
n=β1α1×β2α2×...×βkαk，其中
α
i
>
0
\alpha_i > 0
αi>0。
则
φ
(
n
)
=
n
(
1
−
1
p
1
)
(
1
−
1
p
2
)
.
.
.
(
1
−
1
p
k
)
\varphi (n) = n(1- \frac {1}{p_1})(1-\frac 1 p_2)...(1-\frac 1 p_k)
φ(n)=n(1−p11)(1−p12)...(1−p1k)
（3）证明
所有和n不互质的数，就是n的所有质因子的倍数，剩下的数就一定和n互质，那么问题就转化成了求n的质因子的倍数有多少个 。这里就要用容斥原理来求。我们先用n减去所有质因数的倍数，即
n
−
∑
i
=
1
k
n
p
i
n-\sum_{i = 1} ^k\frac n p_i
n−∑i=1kpni；根据容斥原理，这里多减了所有两个质因数的倍数，因此还要再加上
∑
1
k
n
p
i
p
j
(
i
≠
j
)
\sum_{1} ^k\frac n {p_ip_j} (i\neq j)
∑1kpipjn(i=j)；同时，这里又把同时三个质因数的倍数加回来了，还得减去
∑
1
k
n
p
i
p
j
p
k
(
i
≠
j
≠
k
)
\sum_{1} ^k\frac n {p_ip_jp_k} (i\neq j\neq k)
∑1kpipjpkn(i=j=k) …，即根据容斥原理最后列出的公式为
n
−
∑
i
=
1
k
n
p
i
+
∑
1
k
n
p
i
p
j
−
∑
1
k
n
p
i
p
j
p
k
.
.
.
n-\sum_{i = 1} ^k\frac n p_i +\sum_{1} ^k\frac n {p_ip_j} - \sum_{1} ^k\frac n {p_ip_jp_k}...
n−∑i=1kpni+∑1kpipjn−∑1kpipjpkn...，将原公式
φ
(
n
)
\varphi (n)
φ(n)展开正好就是该公式。
该公式的时间复杂度卡在分解质因数上，是 O n O\sqrt n On
（4）代码
int euler(int x)
{
int ans = x;
for(int i = 2;i <= x/i;i++)
{
if(x % i == 0)
{
ans -= ans/i;
while(x % i == 0)
x /= i;
}
}
if(x > 1)
ans -= ans/x;
return ans;
}
2、筛法求欧拉函数之和
我们已知一个数的欧拉函数公式求解时间复杂度是
O
(
n
)
O(\sqrt n)
O(n)，那么我们求1~n的所有数的欧拉函数之和的时间复杂度就是
O
(
n
n
)
O(n\sqrt n)
O(nn)。然而当我们利用欧拉筛来求n个数的欧拉函数之和能将复杂度降到
O
(
n
)
O(n)
O(n)。
我们记phi[i]为i的欧拉函数。
（1）原理
在欧拉筛中，我们将所有合数标记了一次，所有的质数存了一次，那么在每次标记合数或存质数的时候，我们可以求出他的欧拉函数。
首先，质数的欧拉函数最好求，质数x的欧拉函数就等于x-1。
对于每个合数，我们是标记的primes[j] * i，我们将该数分成两种情况来讨论：
①i % primes[j] == 0，说明primes[j]是i的一个质因数，那么primes[j] * i的质因数和i的质因数完全相同，那么他们的欧拉函数公式仅前面的n不同，因此这种情况的phi[primes[j] * i] = primes[j] * phi[i]。
②i % primes[j] != 0，primes[j] * i的质因数只比i少了primes[j]这一项，其他相同，那么根据公式可得phi[primes[j] * i] = primes[j] * (1-1/primes[j]) * phi[i]化简即phi[primes[j] * i] = (primes[j] - 1) * phi[i]。
（2)代码
int cnt,primes[MAXN],phi[MAXN];
bool vis[MAXN];
int eulers(int n)
{
phi[1] = 1;
for(int i = 2;i <= n/i;i++)
{
if(!vis[i])
{
primes[cnt++] = i;
phi[i] = i-1;
}
for(int j = 0;primes[j] <= n/i;j++)
{
vis[primes[j] * i] = 1;
if(i % primes[j] == 0)
{
phi[primes[j] * i] = primes[j] * phi[i];
break;
}
phi[primes[j] * i] = (primes[j]-1) * phi[i];
}
}
}
2、欧拉定理
若a和n互质，则 a φ ( n ) ≡ 1 a^{\varphi (n)} \equiv 1 aφ(n)≡1 (mod n)
3、费马小定理
费马小定理就是欧拉定理的一个特殊情况
已知a和p互质，且p为质数，则
a
p
−
1
≡
1
a^{p-1}\equiv1
ap−1≡1 (mod p)。
传送门：欧拉定理与费马小定理的证明
四、快速幂
1、定义
若计算a的n次幂是一般要循环n次，让n个a相乘，时间复杂度为 O ( n ) O(n) O(n)，若n很大时如n = 1 0 9 10^9 109 时要循环 1 0 9 10^9 109次。为了优化时间复杂度，我们用快速幂达到 O ( log n ) O(\log n) O(logn)的复杂度。
2、原理
我们将幂次n按二进制分解，分解为
1
、
2
、
4
、
8
、
.
.
.
、
k
1、2、4、8、...、k
1、2、4、8、...、k，使得
1
+
2
+
4
+
.
.
.
+
k
=
n
1+2+4+...+k = n
1+2+4+...+k=n，然后我们直学计算底数
a
1
×
a
2
×
a
4
×
.
.
.
×
a
k
a^1\times a^2\times a^4 \times ... \times a^k
a1×a2×a4×...×ak得到的答案就是
a
n
a^n
an。
这里对幂次的分解用到了位运算，用二进制的特性就很好计算了，详情见此次。
举个例子：如计算37时，我们直接计算
3
1
×
3
2
×
3
4
3^1\times 3^2\times3^4
31×32×34，这样只需要循环三次即可。
3、代码模板：
根据题目中间还可以加上取余运算，如下，设余数是mod。
int qmi(int a,int n) //计算a的n次幂
{
long long res = 1;
while(n)
{
if(n & 1)
res = (1ll)*res * a % mod; //用long long防止乘法溢出
n >>= 1;
a = (1ll)*a * a % mod;
}
return res;
}
五、欧几里得
1、欧几里得算法
（1）定义
欧几里得算法，又称辗转相除法，用于求两个非负整数a、b的最大公约数。
（2）公式
gcd(a,b) = gcd(b,a mod b)
（3）证明
①假设我们要求
g
c
d
(
a
,
b
)
gcd(a,b)
gcd(a,b)，a可以表示成
a
=
k
b
+
r
a = kb + r
a=kb+r（k、r均为正整数），因此a mod b = r。又设d为(a,b)的任意一个公约数。
②
r
=
a
−
k
b
r = a-kb
r=a−kb，两边同时除以d得
r
d
=
a
d
−
k
b
d
\frac r d = \frac a d - \frac {kb} d
dr=da−dkb，a、kb均是d的倍数，因此
r
d
\frac rd
dr是一个整数，即d是余数r的约数。
③而 d 又是 b 的约数，所以d也是（b，a mod b）的公约数。从而证得 （a，b） 的任意公约数也是（b，a mod b）的公约数。
④因此（a，b） 的最大公约数也是（b，a mod b）的最大公约数。证毕。
（4）代码模板
int gcd(int a,int b)
{
if(b)
return gcd(b,a % b); //递归辗转相除
return a; //a和0的最大公约数是a
}
2、拓展欧几里得算法
在看拓展欧几里得算法之前一定要了解裴蜀定理。
（1）裴蜀定理
①定义
对于任意正整数a、b，一定存在整数x和y，使得
a
x
+
b
y
=
g
c
d
(
a
,
b
)
ax+by=gcd(a,b)
ax+by=gcd(a,b)。
②证明
假设
g
c
d
(
a
,
b
)
=
d
gcd(a,b) = d
gcd(a,b)=d，任取整数x、y（也有可能为负数），
a
x
+
b
y
=
k
d
ax+by=kd
ax+by=kd，因为 a 和 b 都是 d 的整数倍，因此
a
x
+
b
y
ax+by
ax+by也肯定是d的整数倍。所以
a
x
+
b
y
ax+by
ax+by能凑出来的最小正整数就是d（最小也是1倍）。
③重要推论
裴蜀定理的一个重要推论是：a和b互质的一个充分必要条件是存在整数x和y使得
a
x
+
b
y
=
1
ax+by = 1
ax+by=1。
（2）拓展欧几里得算法
在了解裴蜀定理之后，那么如何构造一组x、y使得 a x + b y = g c d ( a , b ) ax+by=gcd(a,b) ax+by=gcd(a,b)呢，这就是拓展欧几里得算法的内容。
拓展欧几里得的求法就是在欧几里得算法求gcd的同时，求出a、b的一组x、y。那么就得来看欧几里得算法了，这里求x、y时，参数一定要传地址，因为这里要递归求x和y，所以要传引用。
①首先，当b=0时，此时a、b的gcd等于a，因此此时的x=1，y=0就是一组合法解
a
×
1
+
b
×
0
=
a
a\times 1+b\times 0 = a
a×1+b×0=a。
②b不等于0时，我们会递归求b和a mod b的gcd，他们的系数会发生变化，我们将a mod b可以写作
a
−
⌊
a
b
⌋
×
b
a-\lfloor \frac ab\rfloor\times b
a−⌊ba⌋×b，而
g
c
d
(
a
,
b
)
=
g
c
d
(
b
,
a
%
b
)
gcd(a,b)=gcd(b,a\%b)
gcd(a,b)=gcd(b,a%b)那么递归后系数表达式也可以化简为：
b
x
′
+
(
a
%
b
)
y
′
=
g
c
d
(
b
,
a
%
b
)
bx'+(a\%b)y'=gcd(b,a\%b)
bx′+(a%b)y′=gcd(b,a%b)
b
x
′
+
(
a
−
⌊
a
b
⌋
×
b
)
y
′
=
g
c
d
(
b
,
a
%
b
)
bx'+(a-\lfloor \frac ab\rfloor\times b)y'=gcd(b,a\%b)
bx′+(a−⌊ba⌋×b)y′=gcd(b,a%b)整理a和b的系数可得
a
y
′
+
b
(
x
′
−
⌊
a
b
⌋
×
y
′
)
=
g
c
d
(
b
,
a
%
b
)
ay'+b(x'-\lfloor \frac ab\rfloor\times y')=gcd(b,a\%b)
ay′+b(x′−⌊ba⌋×y′)=gcd(b,a%b)由此可得
x
=
y
′
,
y
=
x
′
−
⌊
a
b
⌋
×
y
′
x=y',y=x'-\lfloor \frac ab\rfloor\times y'
x=y′,y=x′−⌊ba⌋×y′因此我们可以先递归求出
x
′
x'
x′和
y
′
y'
y′后再求得x和y。
为了方便代码，我们一般在递归时会交换
x
′
x'
x′和
y
′
y'
y′，因此计算时x不变，y减去
⌊
a
b
⌋
×
x
\lfloor \frac ab\rfloor\times x
⌊ba⌋×x即可。
代码模板：
int exgcd(int a,int b,int &x,int &y)
{
if(b)
{
int d = exgcd(b,a%b,y,x);
y -= a/b * x;
return d;
}
x = 1,y = 0;
return a;
}
（3）线性同余方程：拓展欧几里得算法的应用
已知a、b、m，求解关于x的线性同余方程： a × x ≡ b a\times x\equiv b a×x≡b (mod m)。
对于该方程，首先进行变形，得：
(
a
×
x
)
%
m
=
b
(a\times x)\%m=b
(a×x)%m=b，那我们设一个整数
y
′
=
⌊
a
×
x
m
⌋
y'=\lfloor \frac {a\times x}m\rfloor
y′=⌊ma×x⌋，可得
a
×
x
−
m
×
y
′
=
b
a\times x-m\times y'=b
a×x−m×y′=b，设
y
=
−
y
′
y=-y'
y=−y′，则原式可化简为
a
x
+
m
y
=
b
ax+my=b
ax+my=b，这不就是拓展欧几里得的形式嘛。
因此，当b是
g
c
d
(
a
,
m
)
gcd(a,m)
gcd(a,m)的倍数时，原式有解，由拓展欧几里得算法可解得a的系数
x
′
x'
x′，则
x
=
b
g
c
d
(
a
,
m
)
×
x
′
x=\frac b{gcd(a,m)}\times x'
x=gcd(a,m)b×x′；否则无解。
六、逆元
1、定义
我们一般要用到的逆元是一个数模m的乘法逆元，一个数a 有逆元的充分必要条件是a必须和m互质。设某数a模m的乘法逆元是x，则 a × x ≡ 1 a\times x\equiv1 a×x≡1 (mod m)。利用乘法逆元我们可以将除法转换成乘法，这样就可以有效的防止爆精度或小数无法取模的情况了，因此逆元的用处很多。
2、逆元的求法
求一个数a模m的乘法逆元时，有两种情况：
①当m为质数时，用快速幂求逆元。a的逆元是
a
m
−
2
a^{m-2}
am−2
②当m不是质数时，用拓展欧几里得求逆元，exgcd(a,p,x,y)求得一个x，然后令
x
=
(
x
+
m
)
%
m
x=(x+m)\%m
x=(x+m)%m，保证了得到的x是一个正整数，则该x就是a的逆元。
3、证明
①若m是质数，根据定义有
n
a
≡
n
×
x
\frac na\equiv n\times x
an≡n×x (mod m)，两边同乘a可得
n
≡
n
×
a
×
x
n\equiv n\times a\times x
n≡n×a×x (mod m)，对于质数m来说一定会存在一个和他互质的数n，因此可以约去两边的n，得
a
×
x
≡
1
a\times x\equiv 1
a×x≡1 (mod m)。
由费马小定理可知，当m为质数时
a
m
−
1
≡
1
a^{m-1}\equiv 1
am−1≡1 (mod m)，拆出一个a可得
a
×
a
n
−
2
≡
1
a\times a^{n-2} \equiv 1
a×an−2≡1 (mod m)。因此可得：当m为质数时，a的逆元等于
a
m
−
2
a^{m-2}
am−2。
②若m不是质数，a有逆元的充要条件是a与p互质，因此
g
c
d
(
a
,
p
)
=
1
gcd(a,p)=1
gcd(a,p)=1，根据定义假设a得逆元是x，则有
a
×
x
≡
1
a\times x\equiv 1
a×x≡1 (mod p)，等价于
a
x
+
p
y
=
1
ax+py=1
ax+py=1，用拓展欧几里得求得一个正整数x即为a的逆元。
4、代码模板
①若m是质数
int qmi(int a,int b,int m)
{
int res = 1;
while(b)
{
res = (1ll) * res * a % m;
a = (1ll) * a * a % m;
b >>= 1;
}
return res;
}
int main()
{
cin >> a >> m;
if(a % m) //因为m是一个质数，因此只要a不是m的倍数，a就和m互质
cout << "no" << endl;
else
cout << qmi(a,m-2,m) << endl;
return 0;
}
②若m不是质数
int exgcd(int a,int b,int &x,int &y)
{
if(b)
{
int d = exgcd(b,a%b,y,x);
y -= a/b * x;
return d;
}
x = 1,y = 0;
return a;
}
int main()
{
int a,m,x,y;
cin >> a >> m;
int d = exgcd(a,m,x,y);
if(d == 1)
cout << (x + p) % p << endl;
else
cout << "no" << endl;
return 0;
}
七、高斯消元
高斯消元就是通过利用线性代数中的初等行列变换来达到在 O ( n 3 ) O(n^3) O(n3)的时间复杂度范围内，求得一个n元一次方程的解。
初等行列变换
①给某一行乘一个非零的数
②交换某两行
③将某行乘以一个数k后加到另一行上
通过这三个变换后的行列式，其解不变。
基本思想
就是将所有的系数和等式右边的常数提出来看作一个
n
×
(
n
+
1
)
n\times (n+1)
n×(n+1)的行列式，然后通过初等行列变换使他变成一个上三角行列式，然后从下往上依次解出每一个
x
i
x_i
xi即可。
{
a
11
⋅
x
1
+
a
12
⋅
x
2
+
a
13
⋅
x
3
=
b
1
a
21
⋅
x
1
+
a
22
⋅
x
2
+
a
23
⋅
x
3
=
b
2
a
31
⋅
x
1
+
a
32
⋅
x
2
+
a
33
⋅
x
3
=
b
3
\left\{ \begin{matrix} a_{11}\cdot x_1+a_{12}\cdot x_2+a_{13}\cdot x_3=b_1 \\ a_{21}\cdot x_1+a_{22}\cdot x_2+a_{23}\cdot x_3=b_2 \\ a_{31}\cdot x_1+a_{32}\cdot x_2+a_{33}\cdot x_3=b_3 \end{matrix} \right.
⎩⎨⎧a11⋅x1+a12⋅x2+a13⋅x3=b1a21⋅x1+a22⋅x2+a23⋅x3=b2a31⋅x1+a32⋅x2+a33⋅x3=b3
可得系数矩阵：
a
11
a
12
a
13
b
1
a
21
a
22
a
23
b
2
a
31
a
32
a
33
b
3
\begin{array}{cccc} a_{11} & a_{12} & a_{13} & b_1 \\ a_{21} & a_{22} & a_{23} & b_2 \\ a_{31} & a_{32} & a_{33} & b_3 \\ \end{array}
a11a21a31a12a22a32a13a23a33b1b2b3
代码模板
#include <iostream>
#include <cmath>
#include <algorithm>
using namespace std;
const int N = 110;
const double eps = 1e-6;
int n;
double a[N][N];
int gauss()
{
int c,r; //c是列，r是行
for(c = 0,r = 0;c < n;c++)
{
int t = r;
for(int i = r;i < n;i++) //找到当前列绝对值最大的一行（防止精度爆炸）
if(fabs(a[i][c]) > fabs(a[t][c]))
t = i;
if(fabs(a[t][c]) < eps) //如果当前行都是0的话，直接去找下一行
continue;
for(int i = c;i < n+1;i++) //将绝对值最大的一行和当前整行交换
swap(a[t][i],a[r][i]);
for(int i = n;i >= c;i--) //把当前行的第c个数变成1，倒着除以第一个数即可
a[r][i] /= a[r][c];
for(int i = r+1;i < n;i++) //确定了r行c列的1之后，将之后所有行的第c列都变为0
if(fabs(a[i][c]) > eps)
for(int j = n;j >= c;j--)
a[i][j] -= a[r][j] * a[i][c]; //该行第一个数为a[i][c],相当于将该行减去已经确定了的第r行的a[i][c]倍（初等行列变换3）
r++; //当前行++
}
if(r < n) //如果r<n说明会存在某列全是0，说明倒着推时会出现断层，只能是无解或无穷多的解
{
for(int i = r;i < n;i++)
if(fabs(a[i][n] > eps)) //r行之后等号左边都是0，如果右边不是0的话，会出现0！=0的情况，无解
return 2;
return 1; //无穷多的解
}
for(int i = n-1;i >= 0;i--) //如果有唯一解的话，从下往上倒着推一遍，倒着循环行
for(int j = i+1;j < n;j++)
a[i][n] -= a[j][n] * a[i][j]; //该行的每个系数乘上他对应的解，再用右边的值减去他就是当前行的解
return 0;
}
int main()
{
cin >> n;
for(int i = 0;i < n;i++)
for(int j = 0;j < n+1;j++)
cin >> a[i][j];
int t = gauss();
if(t == 0) //有唯一解，输出
for(int i = 0;i < n;i++)
printf("%.2lf\n",a[i][n]);
else if(t == 1) //无穷多个解
cout << "Infinite group solutions" << endl;
else //无解
cout << "No solution" << endl;
return 0;
}
八、组合数学
1、求组合数
（1）公式法逆元求组合数
根据组合数公式
C
a
b
=
a
!
b
!
⋅
(
a
−
b
)
!
C_{a}^{b}=\frac{a!}{b!\cdot(a-b)!}
Cab=b!⋅(a−b)!a!可以求组合数，但因为数太大一般会对一个数取模，这时我们就需要用逆元来将除法转换成乘法来提高精度，并且能取模，因此有
a
!
b
!
⋅
(
a
−
b
)
!
=
a
!
⋅
[
b
!
]
−
1
⋅
[
(
a
−
b
)
]
−
1
\frac{a!}{b!\cdot(a-b)!}=a!\cdot [b!]^{-1}\cdot [(a-b)]^{-1}
b!⋅(a−b)!a!=a!⋅[b!]−1⋅[(a−b)]−1这里的-1表示逆元的意思。
如果给n组询问，每次询问给两个数a、b，求组合数
C
a
b
C_{a}^{b}
Cab的值（
1
<
n
<
100000
,
1
<
b
≤
a
<
100000
1<n<100000,1<b≤a<100000
1<n<100000,1<b≤a<100000）。那么此时我们就可以用公式法求组合数，因为我们可以用
O
(
n
log
n
)
O(n\log n)
O(nlogn)的复杂度求出所有数的阶乘和阶乘的逆（这里的n是最大的数的大小），就可以每次直接输出公式结果。
代码模板：
#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;
const int N = 100010,mod = 1e9+7;
int n,a,b,fact[N],infact[N]; //fact[i]存i的阶乘,infact[i]存i的阶乘的逆元
int qmi(int a,int b) //快速幂求逆元
{
int res = 1;
while(b)
{
if(b & 1)
res = (1ll) * res * a % mod;
a = (1ll) * a * a % mod;
b >>= 1;
}
return res;
}
int main()
{
fact[0] = infact[0] = 1; //0的阶乘和逆元都初始化为1
for(int i = 1;i < N;i++) //预求出所有数的阶乘和逆元
{
fact[i] = (1ll) * fact[i-1] * i % mod;
infact[i] = (1ll) * infact[i-1] * qmi(i,mod-2) % mod;
}
cin >> n; //n组询问
while (n -- )
{
cin >> a >> b;
cout << (1ll) * fact[a] * infact[a-b] % mod * infact[b] % mod << endl; //输出计算公式
}
return 0;
}
（2）Lucas卢卡斯定理求组合数
当p很小，而a、b很大时，一般用Lucas定理来求解
C
a
b
C_a^b
Cab
m
o
d
mod
mod
p
p
p。
(
1
≤
n
≤
20
,
1
≤
b
≤
a
≤
1
0
18
,
1
≤
p
≤
1
0
5
)
(1≤n≤20,1≤b≤a≤10^{18},1≤p≤10^5)
(1≤n≤20,1≤b≤a≤1018,1≤p≤105)
结论
C
a
b
≡
C
a
p
b
p
⋅
C
a
%
p
b
%
p
C_{a}^{b}\equiv C_{\frac ap}^{\frac bp}\cdot C_{a\%p}^{b\%p}
Cab≡Cpapb⋅Ca%pb%p (mod p)
证明： Lucas定理的证明
代码模板：
#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long LL;
int p,n;
int qmi(int a,int b) //快速幂求逆元
{
LL res = 1;
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) //公式求C(a,b)
{
int res = 1;
for(int i = 1,j = a;i <= b;i++,j--)
{
res = (LL)res * j % p; //分子从a-b+1乘到a
res = (LL)res * qmi(i,p-2) % p; //分母是b的阶乘，转化为乘以阶乘的逆元
}
return res;
}
int lucas(LL a,LL b) //Lucas定理
{
if(a <p && b < p) //如果a和b都小于p，可以直接用公式法求c
return C(a,b);
return (LL)C(a%p,b%p)*lucas(a/p,b/p) % p; //否则就用Lucas定理，用Lucas结论来求解
}
int main()
{
cin >> n;
while (n -- )
{
LL a,b;
cin >> a >> b >> p; //输入a、b和模数p
cout << lucas(a,b) << endl;
}
return 0;
}
2、卡特兰数
（1）定义
卡特兰数是组合数学中一种常出现于各种计数问题中的数列。
（2）描述
①特殊的01序列：给n个0和n个1，将他们排列成一个序列，要求任意前缀序列中的0的个数不少于1的个数的序列有多少个。
我们可以将这个序列转化到网格坐标图中去，从原点开始，0代表向右走一格，1代表向上走一格，那么每一个01序列对应了网格图中的一条从原点到点（n，n）的路径，而每一条路径又对应着唯一的01序列，因此就可以将原问题求序列转化为一个求路径问题，如下图，而要满足前缀中0的个数不少于1的个数也就是x≥y，也就是说所有路径不能过这条红线。
假设n=6，则所有满足要求的序列就等于从原点到（6，6）的所有路径。例一条合法路径如下图蓝线所示。
那么如何计算从原点到（n，n）的路径条数呢。用组合数学即可，我们一共要走2n步，而在这2n步中要选n步向上走，剩下n步向右走，那么排列组合的个数就等于 C 2 n n C_{2n}^n C2nn。然后我们再减去经过红线的路径条数，即为正确答案，这里就是问题的关键了。
所有经过红线的路径，我们把它和红线的交点之后的路径关于红线做对称。也就是说，在第一次到达红线上的点之后的路做个镜像，如下图。
那么这样一来，所有经过红线的路径最终一定都会到达（5，7）点。因此用到（6，6）的路径数减去到（5，7)的路径数即可。即 a n s = C 12 6 − C 12 5 ans=C_{12}^6-C_{12}^{5} ans=C126−C125
②火车进栈出栈问题，另一种描述是用火车的出栈序列个数来说，一共有n列火车，问有多少种进栈出栈的顺序。
我们可以将一辆进栈表示为0，出栈表示为1，所有火车要进一次出一次，因此总序列长度为2n，而且出栈时栈内一定要有火车，就相当于任意位置，前面的0的个数不能少于1的个数，这就转换成了上一种描述方式，典型的求卡特兰数。
（3）结论
a
n
s
=
C
2
n
n
−
C
2
n
n
−
1
ans = C_{2n}^n-C_{2n}^{n-1}
ans=C2nn−C2nn−1
而将所有的C按公式展开后化简上式可得
a
n
s
=
C
2
n
n
÷
(
n
+
1
)
ans = C_{2n}^n\div (n+1)
ans=C2nn÷(n+1)
（4）适用情况
当我们要求的方案数的每一步只有两个状态时（0或1，进栈或出栈，左括号或右括号…），可以考虑转换为卡特兰数来求。
3、容斥原理
传送门：数论之容斥原理 与经典例题
3195
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
