常用代码模板4——数学知识
数论
质数
在大于1的整数中,如果只有1和本身这两个约数,就被称为质数或素数
所有小于等于1的整数既不是质数也不是合数
试除法判定质数
时间复杂度:\(O(\sqrt n)\)
d|n
代表的含义是d能整除n,这里的|
代表整除,d是除数- 一个合数的约数总是成对出现的,如果d | n,那么 (n / d ) | n,因此我们判断一个数是否为质数的时候,只需要判断较小的那一个数能否整除n就行了,即只需枚举d <= (n / d),即d * d <= n 或 d <= sqrt(n)就行了
- 注意
i * i <= n
可能会溢出,如当 n = 2 ^ 31 - 1 时,46340 ^ 2 < n,但是 46341 ^ 2 > n,溢出后 i * i 会从 -2147483648 开始递增,那条件会一直成立,会导致超时
bool is_prime(int x)
{
if (x < 2) return false;
for (int i = 2; i <= x / i; i ++ )
if (x % i == 0)
return false;
return true;
}
试除法分解质因数
时间复杂度:\(O(\sqrt n)\)(最坏)\(O(log n)\)(最好)
质因数(素因数或质因子)在数论里是指能整除给定正整数的质数。如果一个质数是某个数的因数,那么就说这个质数是这个数的质因数
算术基本定理:任何一个大于1的自然数 N,如果N不为质数,那么N可以唯一分解成有限个质数的乘积,质因子如重复可以用指数表示:\(N = P_1 ^ {a_1}*P_2 ^ {a_2}*P_3 ^ {a_3}*……*P_n ^ {a_n}\),这里\(P_1<P_2<...<P_n\)均为质数,其中指数\(a_i\)是正整数。这样的分解称为 N 的标准分解式
1没有质因子
只有一个质因子的正整数为质数。
对比朴素做法,由于一个正整数最多有一个大于根号n的质因数,所以我们只需要枚举到\(\sqrt n\),再判断最后一个即可,不需要直接枚举到\(n\)
![](/i/l/?n=23&i=blog/3051960/202211/3051960-20221130234303333-1940981293.png)
// 当x%i时, x中已不含2 ~ i-1中的质因数, i能整除x, 说明i中也不含2 ~ i-1的质因数, 所以i一定是质数
void divide(int x)
{
for (int i = 2; i <= x / i; i ++ )
if (x % i == 0)
{
int s = 0;
while (x % i == 0) x /= i, s ++ ;
cout << i << ' ' << s << endl;
}
// 一个正整数最多有一个大于根号n的质因数(反证法),如果此时x>1,说明x = Pn * 1,它本身的值就是最后一个质因数
if (x > 1) cout << x << ' ' << 1 << endl;
cout << endl;
}
筛法求1~n的素数
1. 朴素筛法求素数
![](/i/l/?n=23&i=blog/3051960/202211/3051960-20221130234303682-2080454861.png)
1~n 中包含 m 的倍数(从1倍开始)有 n / m 个
质数定理
1~n中有\(\frac{n}{ln\ n}\)个质数.
朴素筛法:从2开始,删去每个数的倍数,这样保证遇到的数都是质数(每个数被它的所有因数筛一遍)时间复杂度\(O(nlogn)\)
埃氏筛法:在朴素筛法的过程中只用质数项去筛(每个数被它的所有质因数筛一遍)时间复杂度:\(O(nloglogn)\)
// 该模板是埃氏筛法
int primes[N], cnt; // primes[]存储所有素数
bool st[N]; // st[x]存储x是否被筛掉
void get_primes(int n)
{
for (int i = 2; i <= n; i ++ )
{
if (st[i]) continue;
primes[cnt ++ ] = i;
for (int j = i + i; j <= n; j += i)
st[j] = true;
}
}
2. 线性筛法求素数
时间复杂度:\(O(n)\)
原理:每个数只被它的最小质因数筛掉,这样每个数就只筛了一遍
分析:我们保证如果 i 没被筛掉,i 一定是质数。筛数过程:i 是倍数,每一次从小到大遍历已得到的质数,直到质数与 i 的乘积大于 n (之后的质数与 i 的乘积都会比 n 大,只是求1 ~ n的质数,筛掉比 n 大的没意义)。
如何保证所有合数都是被它的最小质因数筛掉的?
对于每次外循环,i 是倍数,我们会在这次循环中筛掉以2 ~ x为最小质因数的,是其最小质因数 i 倍的合数,也就是说 primes[j] * i 的最小质因数一定是primes[j]。当 i % primes[j] != 0 时,primes[j] 比 i 的最小质因数还要小,那它一定是 primes[j] * i 的最小质因数;当 i % primes[j] == 0,由于质数是从小到大枚举的,那么primes[j] 一定是 i 的最小质因数,所以也一定 primes[j] * i 的最小质因数。因为我们是先筛数后break,所以 primes[j] * i 此时已经被筛掉了。如果不break,接下来 primes[j] 就不是 primes[j] * i 的最小质因数了,因为 i 的最小质因数比接下来的primes[j]小。
int primes[N], cnt; // primes[]存储所有素数
bool st[N]; // st[x]存储x是否被筛掉
void get_primes(int n)
{
// 质数从2开始,求1~n的质数,枚举到n
for (int i = 2; i <= n; i ++ )
{
if (!st[i]) primes[cnt ++ ] = i;
for (int j = 0; primes[j] <= n / i; j ++ )
{
st[primes[j] * i] = true;
if (i % primes[j] == 0) break;
}
}
}
约数
约数,又称因数。整数a除以整数b(b≠0) 除得的商正好是整数而没有余数,我们就说a能被b整除,或b能整除a。a称为b的倍数,b称为a的约数。
试除法求所有约数
时间复杂度:\(O(\sqrt n)\)
一个合数的约数总是成对出现的,如果d | n,那么 (n / d ) | n,所以只需要找到1~\(\sqrt n\)中的因数 i ,另一半因数直接 n / i 即可得到
vector<int> get_divisors(int x)
{
vector<int> res;
for (int i = 1; i <= x / i; i ++ )
if (x % i == 0)
{
res.push_back(i);
if (i != x / i) res.push_back(x / i);
}
sort(res.begin(), res.end());
return res;
}
约数个数和约数之和
因数个数定理
因数个数定理:
如果\(N = P_1^{c_1} * P_2^{c_2} * ... *P_k^{c_k}\),N的每个约数 a 都可以写成
$ a = P_1^{a_1} * P_2^{a_2} * ... *P_k^{a_k}$的形式,其中 $ 0 \leq a_i \leq c_i\(,所以因数个数 就是\)a_1 - a_k\(的取法数,即\) (c_1 + 1) * (c_2 + 1) * ... * (c_k + 1)$,因为 由算数基本定理,两个数的质因数分解式不同,两个数就不同。
n 中包含 m 的倍数(从1倍开始)有 n / m 个
思考:1 ~ n 中所有数的因数个数之和?
解:因为因数和倍数是一一对应的,所以说相当于求1~ n 中所有数的倍数之和,那么:1 有 n / 1 个倍数,2 有 n / 2 个倍数,……,n 有 n / n 个倍数,所以倍数之和 = 因数之和 = n * (1 + 1 / 2 + 1 / 3 + ... + 1 / n) ~ n * ln(n),平均来看,每个数的因数个数有 log n 个。
int范围内的数,约数个数最多的数大约有1500个约数
约数之和:
将公式(p1^0 + p1^1 + ... + p1^c1) * ... * (pk^0 + pk^1 + ... + pk^ck)用分配律展开,每一项可以写成$ P_1^{a_1} * P_2^{a_2} * ... *P_k^{a_k}$的形式,显然这是N的一个因数,其中 $ 0 \leq a_i \leq c_i$,所以说一共有 (c1 + 1) * (c2 + 1) * ... * (ck + 1)项,因此公式的结果就是所有因数的和。
如果 N = p1^c1 * p2^c2 * ... *pk^ck
约数个数: (c1 + 1) * (c2 + 1) * ... * (ck + 1)
约数之和: (p1^0 + p1^1 + ... + p1^c1) * ... * (pk^0 + pk^1 + ... + pk^ck)
欧几里得算法(辗转相除法)
时间复杂度:\(O(logn)\)
整除的性质
模除定义
定义1:
定义2:
a = nq + r, 则 a % n = r
所以:a % n = a - nq = a - n * [a / n]
库函数:__gcd(a, b) (注意前面是两个_)
计算最小公倍数:a * b / gcd(a, b)
证明欧几里得算法:
-
证明gcd(a, b) = gcd(b, a % b):
左推右:假如d是a, b的公约数,那么d|a,d|b;
a % b = a - b * [a / b] ,令c = [a / b],则a % b = a - b*c
由整除性质推论:d|(a - b % c),即d|(a % b),所以d也是b和a%b的公约数。所以左右两边的最大公约数相同。
右推左:假如d是b和a % b的公约数,那么d|(a % b),d|b;
a % b = a - b*c,a = b*c + (a % b)
由整除性质推论:d|a,所以d也是a和b的公约数。
所以左右两边的最大公约数相同。
-
为什么后一个数最终变为0:
-
gcd(n, 0) = n
int gcd(int a, int b)
{
return b ? gcd(b, a % b) : a;
}
欧拉函数
互质
定义:
性质:
\(\exists x, y \in Z, \text{使得 }xa+yb=1\)
判别方法:
- 两个不同的素数一定互质。例如,2与7、13与19。
- 2.一个素数,另一个不为它的倍数,这两个数互质。例如,3与10、5与26。
- 1和任何一个自然数都互质。如1和9908.
- 相邻两个奇数互质。如49与51。
- 相邻两个自然数互质。如15与16。
- 较大数是素数,则两个数互质。如97与88。
- 两数都是合数(二数差较大),较小数所有的素因数, 都不是较大数的因数,这两个数互质。如357与715, 357=3x7x17, 而3、7和17都不是715的因数,故这两数互质。
- 两数都是合数(二数差较小),这两数之差的所有素因数都不是较小数的因数,这两个数互质。如85和78。 85 - 78=7, 7不是78的因数,故这两数互质。
- 两数都是合数,较大数除以较小数的余数(大于"1” )的所有素因数,都不是较小数的因数,则两数互质。如462与221, 462%221=20, 20=2x2x5。 2、 5都不是221的因数,故这两数互质。
- 两数最大公因数是1。
欧拉函数定义
容斥原理定义
证明:
欧拉函数性质
证明:
性质1:若gcd(n, x) = 1, 则gcd(n, n - x) = 1, 所以与 n 互质的数成对出现,有 phi(n) / 2 对,每对和为 n ,所以结果为 n * phi(n) / 2
性质2:直接用欧拉函数定义展开即可
1∼n 中每个数的欧拉函数之和的结果是n ^ 2级别的,注意数据范围可能结果会爆int
求欧拉函数
时间复杂度:\(O(\sqrt n)\)
int phi(int x)
{
int res = x;
for (int i = 2; i <= x / i; i ++ )
if (x % i == 0)
{
res = res / i * (i - 1); // res = res * (1 - 1 / i) 中的小数会直接得0,结果不对,公式要等价变形一下
while (x % i == 0) x /= i;
}
if (x > 1) res = res / x * (x - 1);
return res;
}
筛法求欧拉函数
用于求 1 ~ N 中每个数的欧拉函数:时间复杂度:\(O(\sqrt n)\)
正确性证明:
int primes[N], cnt; // primes[]存储所有素数
int euler[N]; // 存储每个数的欧拉函数
bool st[N]; // st[x]存储x是否被筛掉
void get_eulers(int n)
{
euler[1] = 1;
for (int i = 2; i <= n; i ++ )
{
if (!st[i])
{
primes[cnt ++ ] = i;
euler[i] = i - 1;
}
for (int j = 0; primes[j] <= n / i; j ++ )
{
int t = primes[j] * i;
st[t] = true;
if (i % primes[j] == 0)
{
euler[t] = euler[i] * primes[j];
break;
}
euler[t] = euler[i] * (primes[j] - 1);
}
}
}
同余与剩余系
简化剩余系的乘法的封闭性:如果 a, b 都属于 m 的简化剩余系(a, b是简化剩余系中的同余类的余数,称它们属于 m 的简化剩余系),则 a * b mod m 也属于 m 的简化剩余系,推广得:一个数的简化剩余系在模 n 的意义下都是相同的
证明:gcd(a, m) = 1, gcd(b, m) = 1 -> gcd(a*b, m) = 1
因为 a * b % m = a * b - m * c,所以 a * b % m 的因数一定是 a * b 和 m * c的公因数。由于 a * b 中不含 m 的因数,所以 a * b % m中也不会有 m 的因数(除 1 以外),所以 a * b % m 与 m 互质。也可以用反证法:如果 ab % m = ab - mc 与 m 不互质,说明 ab - mc 中含有 n 的大于 1 的因数,那么 ab = ab - mc + mc中肯定也含有,与 ab 与 n 互质矛盾,证毕。又因为 0 < a * b % m < m,所以 a * b % m 属于 m 的简化剩余系。
可以推广,如果 n 与 m 互质,则 n % m 也与 m 互质
模算数的运算定律
同余性质
模除分配律
费马小定理
欧拉定理
费马小定理是欧拉定理中 n 为质数的情况 (phi[n] = n - 1)
欧拉定理证明:
快速幂
求a^k mod p,时间复杂度 O(logk)。
数据范围 1 <= a , p <= 2e9,0 <= k <= 2e9
注意:p如果很大,相乘的结果会在%p之前爆int,所以要开long long
正确性证明:
// res 存结果
// t 存a^2^i, 每次都要平方取模, 初始值为a^2^0=a
// 如果k=0, 不会进入while循环, 所以要先把res初始化为1%p,这样, 如果k=0 && p>1, 那结果是1, 正确; 如果k=0 && p=1, 结果是0, 正确; 如果k!=0 && p=1, 则res=0, 循环后结果还是0, 正确; 如果k!=0 && p>1, 结果不影响, 正确
int qmi(int a, int k, int p)
{
int res = 1 % p, t = a;
while (k)
{
if (k & 1) res = res * t % p; // 这里可能爆int
t = t * t % p; // 这里可能爆int
k >>= 1;
}
return res;
}
模逆元
模逆元也称为模倒数,或者模逆元。
b 存在模逆元\(b^{-1}\)的充要条件是 b 与模数 n 互质。
\((ab)^{-1}=a^{-1}\times b^{-1}\)
当模数 n 为质数时,\(b^{m-2}\)即为 b 的模逆元,用快速幂求解
证明:(注:如果 b 不是 p 的倍数,又因为 p 是质数,所以此时 b 和 p 的公因数只有1,所以 b 和 p 互质)(b < p 不是必要的)
// 求b模p的乘法逆元
if (b % p == 0) 无解
else 解 = qmi(b, p - 2, p);
龟速幂
在数据爆long long时,保证求出a * k % p的结果,时间复杂度\(O(logk)\)
原理:
typedef long long LL;
LL slow_mul(LL a, LL k, LL p)
{
bool flag = true;
if (k < 0) k = -k, flag = false;
LL ans = 0;
while (k)
{
if (k & 1) ans = (ans + a) % p;
k >>= 1;
a = (a + a) % p;
}
if (flag) return ans;
else return -ans;
}
扩展欧几里得算法
求x, y,使得ax + by = gcd(a, b)
贝祖定理
注:贝祖定理又名裴蜀定理
gcd(a, b)是a, b的线性组合 sa + tb 能得到的最小正整数
正确性证明:
通解:
\[\text{已知 }ax_0+by_0=(a, b)\\ \text{通解为:} \begin{cases}x=x_0-\frac{b}{d}k,& k\in Z\\ y=y_0+\frac{a}{d}k \end{cases} \]int exgcd(int a, int b, int &x, int &y)
{
if (!b)
{
x = 1; y = 0;
return a;
}
int d = exgcd(b, a % b, y, x);
y -= (a/b) * x;
return d;
}
模除恒等式
线性同余方程
拓展欧几里得算法解线性同余方程证明:
当b数据范围到2e9时,很容易爆int,先转成LL,再模m,结果依然正确
因为根据模除恒等式:ans = (LL)x * b / d
ans % m = ans % m % m
// 给定a, m, b, 求 x, 使a * x = b (mod m)
int d = exgcd(a, m, x, y);
if (b % d == 0) 解 = (LL)x * b / d % m
else 无解
中国剩余定理
定理内容
正确性证明:
模板来自中国剩余定理 - OI Wiki (oi-wiki.org)
LL crt(int n, LL * a, LL * m)
{
LL M = 1, ans = 0;
for (int i = 1; i <= n; i ++ ) M *= m[i];
for (int i = 1; i <= n; i ++ )
{
LL Mi = M / m[i], x, y;
exgcd(Mi, m[i], x, y);
ans = (ans + a[i] * Mi * x % M) % M;
}
return (ans % M + M) % M;
}
拓展中国剩余定理
给定 \(n\) 组非负整数 \(a_i, b_i\) ,求解关于 \(x\) 的方程组的最小非负整数解。
\[\begin{cases} x \equiv a_1\ ({\rm mod}\ m_1) \\ x\equiv a_2\ ({\rm mod}\ m_2) \\ ... \\ x \equiv a_n\ ({\rm mod}\ m_n)\end{cases} \]解法:合并方程
注意:k1 = k1 * (a2 - a1) / d
这一步可能爆long long ,可以用龟速幂算slow_mul(k1, (a2 - a1) / d, m2 / d)
// a[]存余数,m[]存模数
LL excrt(int n, int * a, int * m)
{
LL a1 = a[1], m1 = m[1];
bool flag = true;
for (int i = 2; i <= n; i ++ )
{
LL a2 = a[i], m2 = m[i], k1, k2;
LL d = exgcd(m1, m2, k1, k2);
if ((a2 - a1) % d)
{
flag = false;
break;
}
k1 = k1 * (a2 - a1) / d;
LL t = m2 / d;
k1 = (k1 % t + t) % t;
a1 = k1 * m1 + a1; // 注意这两步不要写反
m1 = m1 / d * m2;
}
if (flag) return (a1 % m1 + m1) % m1;
else return -1;
}
高斯消元
1.解线性方程组
时间复杂度\(O(n^3)\)
思路:
-
目标:将增广矩阵化为最简阶梯型矩阵
-
步骤:
-
从左向右枚举系数矩阵的列。
-
对于每一列,先从未确定方程的行中找到这一列中绝对值最大的数所在的行,并将该行换到未确定方程的行的最上面。
-
将该行的第一个数消成1
-
用该行将未确定方程的行的这一列都消成0
-
迭代n次后就得到阶梯型增广矩阵
-
判断解的情况
-
如果有唯一解,从最后一行开始,将系数矩阵对角线以上的数都消成0(对于第 i 行,a[i][0]~a[i][i - 1] 为0;a[i][i] = 1; a[i][i + 1]~a[i][n - 1]为非零。要消去非零,在消的过程中,解向量就会变成唯一解)
-
问题:为什么要找绝对值最大的那一行?
答:一是因为可以寻找一个非零行,二是作浮点数除法时,显然除以一个大数会让精度更好,误差更小。
解m个方程, n 个未知数的话,循环行时边界用m, 循环列时边界用n
// a[N][N]是增广矩阵
int gauss()
{
int c, r;
// 注意for循环里不要重新定义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) continue;
for (int i = c; i <= n; i ++ ) swap(a[t][i], a[r][i]); // 将绝对值最大的行换到最顶端
for (int i = n; i >= c; i -- ) a[r][i] /= a[r][c]; // 将当前行的首位变成1
for (int i = r + 1; i < n; i ++ ) // 用当前行将下面所有行的这一列消成0
if (fabs(a[i][c]) > eps)
for (int j = n; j >= c; j -- )
a[i][j] -= a[r][j] * a[i][c];
r ++ ;
}
if (r < n) // 系数矩阵对角线上不是全1,有右移,不是上三角矩阵
{
for (int i = r; i < n; i ++ )
if (fabs(a[i][n]) > eps)
return 2; // 无解
return 1; // 有无穷多组解
}
for (int i = n - 1; i >= 0; i -- )
for (int j = i + 1; j < n; j ++ )
a[i][n] -= a[i][j] * a[j][n];
return 0; // 有唯一解
}
2.解异或线性方程组
时间复杂度:\(O(n^3)\)
同解变形:一行异或另一行,解不变
思路:
-
目标:将增广矩阵化为最简阶梯型矩阵
-
步骤:
- 从左向右枚举系数矩阵的列。
- 对于每一列,先从未确定方程的行中找到这一列中的第一个1所在的行,并将该行换到未确定方程的行的最上面。
- 用该行将未确定方程的行的这一列中的1都消成0(两行做异或运算)
- 迭代n次后就得到阶梯型增广矩阵
- 判断解的情况
- 如果有唯一解,从最后一行开始,将系数矩阵对角线以上的数都消成0(对于第 i 行,a[i][0]~a[i][i - 1] 为0;a[i][i] = 1; a[i][i + 1]~a[i][n - 1]为0或1,要消去1,在消的过程中,解向量就会变成唯一解)
例子:
![](/i/l/?n=23&i=blog/3051960/202211/3051960-20221130234304357-499678284.jpg)
int gauss()
{
int c, r;
for (c = 0, r = 0; c < n; c ++ )
{
int t = r;
for (int i = r; i < n; i ++ )
if (a[i][c])
{
t = i;
break;
}
if (!a[t][c]) continue; // 这一列全0,跳过看下一列
// 换到最上面
for (int i = c; i <= n; i ++ )
swap(a[t][i], a[r][i]);
// 把下面所有行的这一列的1都消成0(做异或运算)
for (int i = r + 1; i < n; i ++ )
if (a[i][c])
for (int j = c; j <= n; j ++ )
a[i][j] ^= a[r][j];
r ++ ;
}
if (r < n)
{
for (int i = r; i < n; i ++ )
if (a[i][n]) return 2; // 无解 (0 = 1)
return 1; // 无穷多组解
}
for (int i = n - 1; i >= 0; i -- )
for (int j = i + 1; j < n; j ++ )
if (a[i][j])
a[i][n] ^= a[j][n];
return 0;
}
求组合数
递推法求组合数
时间复杂度\(O(n^2)\),le3级别
公式:\(\binom{a}{b}=\binom{a-1}{b}+\binom{a-1}{b-1}\)
只需把所有数初始化为0,j = 0是结果得1,就可以用公式算出所有组合数
// c[a][b] 表示从a个苹果中选b个的方案数
for (int i = 0; i < N; i ++ )
for (int j = 0; j <= i; j ++ )
if (!j) c[i][j] = 1;
else c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % mod;
通过预处理逆元的方式求组合数
时间复杂度\(O(nlogk)\),1e5级别,a , b < p(如果a , b > p,能会出现a! % p = 0的情况,而在a!中的p可能在组合数公式中是会被分母消掉的)
原理:
- 模数固定,多组a, b
首先预处理出所有阶乘取模的余数fact[N],以及所有阶乘取模的逆元infact[N]
如果取模的数是质数,可以用费马小定理求逆元
int qmi(int a, int k, int p) // 快速幂模板
{
int res = 1;
while (k)
{
if (k & 1) res = (LL)res * a % p;
a = (LL)a * a % p;
k >>= 1;
}
return res;
}
// 预处理阶乘的余数和阶乘逆元的余数
fact[0] = infact[0] = 1;
for (int i = 1; i < N; i ++ )
{
fact[i] = (LL)fact[i - 1] * i % mod;
infact[i] = (LL)infact[i - 1] * qmi(i, mod - 2, mod) % mod;
}
- 模数不固定
a, b为1e5级别, a, b < p, p为质数, 求C(a, b) % p
int C(int a, int b, int p) // 通过定理求组合数C(a, b)
{
if (a < b) return 0;
LL x = 1, y = 1; // x是分子,y是分母
for (int i = a, j = 1; j <= b; i --, j ++ )
{
x = (LL)x * i % p;
y = (LL) y * j % p;
}
return x * (LL)qmi(y, p - 2, p) % p;
}
Lucas定理求组合数
时间复杂度\(O(p\log_{2}n)\),p是1e5级别,n是1e18级别, a , b > p
lucas定理:
若p是质数,则对于任意整数 1 <= m <= n,有:
C(n, m) = C(n % p, m % p) * C(n / p, m / p) (mod p)
证明:
C函数原理:
int qmi(int a, int k, int p) // 快速幂模板
{
int res = 1 % p;
while (k)
{
if (k & 1) res = (LL)res * a % p;
a = (LL)a * a % p;
k >>= 1;
}
return res;
}
int C(int a, int b, int p) // 通过定理求组合数C(a, b)
{
if (a < b) return 0;
LL x = 1, y = 1; // x是分子,y是分母
for (int i = a, j = 1; j <= b; i --, j ++ )
{
x = (LL)x * i % p;
y = (LL) y * j % p;
}
return x * (LL)qmi(y, p - 2, p) % p;
}
int lucas(LL a, LL b, int p)
{
if (a < p && b < p) return C(a, b, p);
return (LL)C(a % p, b % p, p) * lucas(a / p, b / p, p) % p;
}
分解质因数法求组合数
当我们需要求出组合数的真实值,而非对某个数的余数时,分解质因数的方式比较好用:
1. 筛法求出范围内的所有质数
2. 通过 C(a, b) = a! / b! / (a - b)! 这个公式求出每个质因子的次数。 n! 中p的次数是 n / p + n / p^2 + n / p^3 + ...
3. 用高精度乘法将所有质因子相乘
int primes[N], cnt; // 存储所有质数
int sum[N]; // 存储每个质数的次数
bool st[N]; // 存储每个数是否已被筛掉
void get_primes(int n) // 线性筛法求素数
{
for (int i = 2; i <= n; i ++ )
{
if (!st[i]) primes[cnt ++ ] = i;
for (int j = 0; primes[j] <= n / i; j ++ )
{
st[primes[j] * i] = true;
if (i % primes[j] == 0) break;
}
}
}
int get(int n, int p) // 求n!中的次数
{
int res = 0;
while (n)
{
res += n / p;
n /= p;
}
return res;
}
vector<int> mul(vector<int> a, int b) // 高精度乘低精度模板
{
vector<int> c;
int t = 0;
for (int i = 0; i < a.size(); i ++ )
{
t += a[i] * b;
c.push_back(t % 10);
t /= 10;
}
while (t)
{
c.push_back(t % 10);
t /= 10;
}
return c;
}
get_primes(a); // 预处理范围内的所有质数
for (int i = 0; i < cnt; i ++ ) // 求每个质因数的次数
{
int p = primes[i];
sum[i] = get(a, p) - get(b, p) - get(a - b, p);
}
vector<int> res;
res.push_back(1);
for (int i = 0; i < cnt; i ++ ) // 用高精度乘法将所有质因子相乘
for (int j = 0; j < sum[i]; j ++ )
res = mul(res, primes[i]);
卡特兰数
该问题可以转换成从坐标(0,0)走到(n,n)的合法路径数
因为:
- 序列每一位只有0和1两种选择,坐标图中每一步也是向右向上两种走法。
- 序列中0和1都要求n个,坐标图中做到(n,n)也要求向右向上都是n步,所以两者方法总数都是C(2n, n)
转换之后分析:合法(x >= y),则若路径经过图中红线就不合法。又因为经过红线的路径经过轴对称可以证明一定会到达(n - 1, n + 1), 且由(0,0)走到(n - 1, n + 1)一定会经过红线,所以两者一一对应。所以合法路径数:
C(2n, n) - C(2n, n - 1) = C(2n, n) / (n - 1)
给定n个0和n个1,它们按照某种顺序排成长度为2n的序列,满足任意前缀中0的个数都不少于1的个数的序列的数量为:
Cat(n) = C(2n, n) / (n + 1)
容斥原理
时间复杂度\(O(2^n)\)
由二项式定理,一共 \(2^n-1\) 项
证明:C(n, 1) + C(n, 2) + ... + C(n, n) = (1 + 1)^n - C(n, 0) = 2^n - 1
表示s[1 ~ n]中选择的状态的方法:
-
dfs,搜索树有 n 层,每一层有两种情况,是否选择s[i]
-
用数的二进制表示。一个数的二进制串中,从第 0 位开始,第 i 位为 1 ,表示选择了s[i],第 i 位为 0,表示不选择s[i]。
状态表示的数为[1 ~ (1 << n) - 1] (从第0位到第n-1位全1的十进制大小为2^n - 1,1 << n == 2 ^ n)
简单博弈论
NIM游戏
给定N堆物品,第i堆物品有Ai个。两名玩家轮流行动,每次可以任选一堆,取走任意多个物品,可把一堆取光,但不能不取。取走最后一件物品者获胜。
我们把这种游戏称为NIM博弈。把游戏过程中面临的状态称为局面。整局游戏第一个行动的称为先手,第二个行动的称为后手。
每一个局面只有两种状态:
-
(先手)必败状态:在当前局面,先行动者会输掉游戏。即不存在某种方式,将当前局面变成先手必败状态,也就是说,无论怎么行动,局面都会变成先手必胜状态,那对方必胜、我方必败。
-
(先手)必胜状态:存在某种方式,将当前局面变成先手必败状态。由于之后的先手必败局面是对手先行动,所以对手必败、我方必胜。
所谓采取最优策略是指,若在某一局面下存在某种行动,使得行动后局面变成先手必败局面,则优先采取该行动。
我们讨论的博弈问题一般都只考虑理想情况,即两人均无失误,都采取最优策略行动时游戏的结果,因此先手必胜和先手必败会交替出现,直到游戏结束。
异或运算律
定理: NIM博弈先手必胜,当且仅当 A1 ^ A2 ^ ... ^ An != 0
证明:
依据证明
0 ^ 0 ^ ... ^ 0 = 0,一定是一个先手必败态。
如果当前A1 ^ A2 ^ A3 ^ ... ^ An != 0,则一定存在某种方式,使得A1 ^ A2 ^ A3 ^ ... ^ An == 0
证明:设异或的结果x的二进制表示中,1的最高位是第k位,那么在A1到An中,一定存在Ai,Ai的第k位为1。
这里可以使用反证法,如果A1到An中第k位都是0,则异或的结果x的第k位一定是0,矛盾,所以Ai一定存在。
Ai ^ x 一定小于Ai。Ai ^ x的第k位一定是0,那么即使Ai ^ x的其余位都是1,结果也是2 ^ k -1, 小于2 ^ k。而Ai >= 2 ^ k
所以如果令Ai减少到Ai ^ x,那异或的结果A1 ^ A2 ^ ... ^ Ai ^ x ^ ... ^ An == x ^ x == 0。证毕。如果当前A1 ^ A2 ^ A3 ^ ... ^ An == 0,则无论怎么做,一定会得到A1 ^ A2 ^ A3 ^ ... ^ An != 0的结果。
证明:利用反证法:如果令Ai -> Aj,使得A1 ^ A2 ^ ... ^ Aj ^ x ^ ... ^ An == 0。该式与变化前的A1 ^ A2 ^ ... ^ Ai ^ x ^ ... ^ An == 0异或,由异或的归零律得Ai ^ Aj == 0,则Ai == Aj,即Ai不变,矛盾,所以原命题成立。
结合问题
如果当前局面异或结果不为0,则先手一定可以令异或结果为0,则后手再操作后一定会使异或结果为0。依次类推可得,如果当前局面异或结果不为0,则先手可以一直不为0,后手一直为0,直到所有的A都变为0,后手必败,先手必胜;如果当前局面异或的结果为0,则先手一定会使异或结果不为0,后手一定可以是异或结果变为0。依次类推,先手一直是0,后手一直不为0,最终先手面对全0结果,先手必败,后手必胜。
公平组合游戏ICG
若一个游戏满足:
-
由两名玩家交替行动;
-
在游戏进程的任意时刻,可以执行的合法行动与轮到哪名玩家无关;
-
不能行动的玩家判负;
则称该游戏为一个公平组合游戏。
NIM博弈属于公平组合游戏,但城建的棋类游戏,比如围棋,就不是公平组合游戏。因为围棋交战双方分别只能落黑子和白子,胜负判定也比较复杂,不满足条件2和条件3。
有向图游戏
给定一个有向无环图,图中有一个唯一的起点,在起点上放有一枚棋子。两名玩家交替地把这枚棋子沿有向边进行移动,每次可以移动一步,无法移动者判负。该游戏被称为有向图游戏。
任何一个公平组合游戏都可以转化为有向图游戏。具体方法是,把每个局面看成图中的一个节点,并且从每个局面向沿着合法行动能够到达的下一个局面连有向边。
必胜状态:可以走到某一个必败状态
必败状态:走不到任何一个必败状态 / 只能走到必胜状态
Mex运算
设S表示一个非负整数集合。定义mex(S)为求出不属于集合S的最小非负整数的运算,即:
mex(S) = min{x}, x属于自然数,且x不属于S
也就是能到达的最小自然数。
SG函数
在有向图游戏中,对于每个节点x,设从x出发共有k条有向边,分别到达节点y1, y2, …, yk,定义SG(x)为x的后继节点y1, y2, …, yk 的SG函数值构成的集合再执行mex(S)运算的结果,即:
SG(x) = mex({SG(y1), SG(y2), …, SG(yk)})
特别地,整个有向图游戏G的SG函数值被定义为有向图游戏起点s的SG函数值,即SG(G) = SG(s),终点的SG值定义为0。
有向图游戏及其和的SG分析
设G1, G2, …, Gm 是m个有向图游戏。定义有向图游戏G,它的行动规则是任选某个有向图游戏Gi,并在Gi上行动一步。G被称为有向图游戏G1, G2, …, Gm的和。
有向图游戏的和的SG函数值等于它包含的各个子游戏SG函数值的异或和,即:
SG(G1, G2, ..., Gm) = SG(G1) ^ SG(G2) ^ … ^ SG(Gm)
定理:
单个有向图游戏
![]()
有向图游戏的某个局面必胜, 当且仅当该局面对应节点的SG函数值大于0。
有向图游戏的某个局面必败, 当且仅当该局面对应节点的SG函数值等于0。
证明:
终点的SG值等于0,为先手必败态。
如果该局面的SG函数值大于0,则它一定可以到达某个SG函数值等于0的局面。
因为该节点的SG值是不在子节点SG值中的最小自然数,如果没有SG值等于0的子节点,则该节点的SG值应该取0,矛盾,因此一定有SG值等于0的子节点。
如果该局面的SG函数值等于0,则它一定不能到达SG函数值等于0的局面。
因为如果2已经证明,如果存在SG值等于0的子节点,则该节点的SG值一定大于0
综上,结合有向图游戏局面的两种状态定义。SG值等于0为先手必败态,SG值大于0为先手必胜态。
有向图游戏的和
有向图游戏的和的某个局面必胜,当且仅当所有子游戏的当前局面的SG函数值的异或和大于0。(一般直接取起点的SG函数值异或和判断)
有向图游戏的和的某个局面必败,当且仅当所有子游戏的当前局面的SG函数值的异或和等于0。(一般直接取起点的SG函数值异或和判断)
证明:
所有子游戏的终点的SG值都为0,0 ^ 0 ^ ... ^ 0 = 0,为先手必败态。
如果当前有向图游戏的和的异或结果大于0,则一定存在某种方式,使得异或和变成等于0。
证明:设异或的结果x的二进制表示中,1的最高位是第k位,那么在SG1到SGn中,一定存在SGi,SGi的第k位为1。
这里可以使用反证法,如果SG1到SGn中第k位都是0,则异或的结果x的第k位一定是0,矛盾,所以SGi一定存在。
SGi ^ x 一定小于SGi。SGi ^ x的第k位一定是0,那么即使SGi ^ x的其余位都是1,结果也是2 ^ k -1, 小于2 ^ k。而SGi >= 2 ^ k
由于SGi是不在子节点SG值中的最小自然数,那么一定存在SG值为 0,1, ... , SGi - 1的子节点(反证法),所以一定可以令SGi减小到SGi ^ x。那异或的结果就会变成0。证毕。如果当前有向图游戏的和的异或结果等于0,无论怎么走,一定会使异或结果变成大于0。(反证法)
结合NIM游戏的推导,得证。
标签:...,int,代码,异或,数学知识,primes,SG,质数,模板 From: https://www.cnblogs.com/zhengzirui/p/16940194.html