数学填坑专辑
裴蜀定理
拓展欧几里得定理
当已知\(a\)，\(b\)时，求一组\(x\)，\(y\)使得\(a\times x + b\times y = GCD(a, b)\)
-
解一定存在
-
因为\(GCD(a,b) = GCD(b, a \% b)\)，所以\(a \times x + b \times y = GCD(b, a \%b ) = b \times x + (a \% b) \times y = b \times x + (a - \lfloor \frac ab \rfloor \times b) \times y = a \times y + b \times (x - \lfloor \frac ab \rfloor \times y)\)。
所以，我们成功地由\(a\)与\(b\)的线性组合转化成了\(b\)与\(a \%b\)的线性组合。
说人话，我们可以在求\(GCD(a,b)\)的过程中求出一组\(x\)和\(y\)了。
应用：
-
求解线性同余方程
定理1：对于方程\(a \times x + b \times y = c\)（等价于\(a \times x \equiv c\pmod b\)）有整数解的充要条件为，\(GCD(a,b)\mid c\)。（裴蜀定理）
根据定理一，我们可以用拓展欧几里得定理求出\(a \times x + b \times y = GCD(a,b)\)的一组解，再\(\div GCD(a,b)\)并\(\times c\)就行了。
-
求上面那个方程的最小解
定理2：\(x_i = x_0 + \frac{b}{GCD(a,b)} \times t\)。其中，\(x_i\)为不定方程的任意解，\(x_0\)为已知解，\(t\)为任意正整数。
由定理2可知，\(x_{min} = x_0 + \frac{b}{GCD(a,b)} \times t\)，所以\(x_{min} = x_0 \mod \frac{b}{GCD(a,b)}\)。别忘了特判\(x \leq 0\)的情况，若\(x \leq 0\)，\(x_{min} = x_0 \mod \frac{b}{GCD(a,b)} + \frac{b}{GCD(a,b)}\)。
（其实上面那个式子是看花姐的，本人太弱了，并不会证……
整数分解
题意：给出两个整数\(a\)和\(b\)，求\(a^b\)的因子和。
前置芝士：整数的唯一分解定理，约数和公式，等比数列求和，快速幂
-
整数的唯一分解定理
任意正整数都有且只有一种方式写出其素因子的乘积表达式：
\[A = (p_1 ^ {k_1}) \times (p_{2} ^ {k_2}) \times (p_3 ^{k_3} * ) \times \cdot\cdot\cdot \times(p_n^{k_n}) \]其中，\(p_i\)为素数。
-
约数和公式
对于已经分解的整数\(A = (p_1 ^ {k_1}) \times (p_{2} ^ {k_2}) \times (p_3 ^{k_3} * ) \times \cdot\cdot\cdot \times(p_n^{k_n})\)
\(A\)的所有因子之和为：
\[S = (1+p_1+p_1^2+p_1^3+\cdot\cdot\cdot+p_1^{k_1})\times(1+p_2+p_2^2+p_2^3+\cdot\cdot\cdot+p_2^{k_2})\times\cdot\cdot\cdot\times(1+p_n+p_n^2+p_n^3+\cdot\cdot\cdot+p_n^{k_n}) \] -
等比数列求和(用递归二分求\(1+p+p^2+p^3+\cdot\cdot\cdot+p^n\))
-
若\(n\)为\(0\)，则答案为\(1\)。
-
若\(n\)为奇数，一共有偶数项，则
\[1+p+p^1+p^2+p^3+\cdot\cdot\cdot+p^n\\=(1+p^{n/2+1})+p\times(1+p^{n/2+1})+\cdot\cdot\cdot+p^{n/2}\times(1+p^{n/2+1})\\=(1+p+p^2+\cdot\cdot\cdot+p^{n/2})\times(1+p^{n/2+1}). \] -
若\(n\)为偶数，就可以把\(n\)看做\(n-1\)，按照奇数做法，将最后答案\(+p^n\)就可以了
-
-
快速幂（过于基础，不再赘述）
然后我们就可以愉快地做这道题了：）
参考代码
#include <cstdio>
#define LL long long
using namespace std;
const int mode = 9901, maxn = 1e4 + 10;;
LL n[maxn],p[maxn],cnt,a,b;
LL q_pow(LL x, int y){
LL ans = 1;
while(y){
if(y & 1) ans *= x, ans %= mode;
x *= x, x %= mode;
y >>= 1;
} return ans;
}
int get_sum(LL p, LL n){
if(n == 0) return 1;
LL ans = 0;
if(!(n & 1)) ans = q_pow(p, n), n -= 1;
return (ans + (get_sum(p, n / 2) * (1 + q_pow(p, n / 2 + 1)) % mode) % mode ) % mode;
}
int main(){
scanf("%lld%lld", &a, &b);
for(int i = 2; i * i <= a;){ //根号法+递归法
if(a % i == 0){
p[++cnt] = i;
n[cnt] = 0;
while(a % i == 0) n[cnt]++, a /= i;
}
if(i == 2) i += 1; //奇偶法
else i += 2;
}
if(a != 1) p[++cnt] = a, n[cnt] = 1;
LL Ans = 1;
for(int i = 1; i <= cnt; ++ i)
Ans = (Ans * get_sum(p[i], n[i] * b)) % mode; // 别忘了将次数*b
printf("%lld\n", Ans);
return 0;
}
做法源于数学一本通（不然我也想不出这么奇怪的命名方式……）
如何\(n\log n\)分解\(n!\):
先筛出\([1,n]\)中每一个质数\(p\)，然后考虑\(n!\)中一共有多少个质因子\(p\)，\(n!\)中质因子\(p\)的个数等于\([1,n]\)中每一个数质因子中\(p\)的个数之和。在\([1,n]\)中，至少包含一个质因子\(p\)的有\(n/p\)个，至少包含两个质因子\(p\)的有\(n/p^2\)个，但是其中的一个因子在包含一个的时候算过了，所以只需要另外统计一个，即累加上\(n/p^2\)。总复杂度\(O(n\log n)\)。
组合数计算
-
加法递推：\(C(n, m) = C(n - 1, m) + C(n - 1, m - 1)\)
复杂度：\(O(n^2)\)
-
乘法递推：\(C(n,m)=C(n,m-1)*(n-m+1)/m\)
复杂度：\(O(n)\)（注意要先乘后加！）
-
当模数不大时，可以使用卢卡斯定理来降低复杂度
代码：
int Lucas(int n, int m){ if(!m) return 1; return (C(n % P, m % P) + Lucas(n / P, m / P)) % P; }
以上算法都很容易爆long long，必要时请使用高精度运算。
第二类斯特林数
-
定义：第二类斯特林数\(S(n,m)\)表示的是把\(n\)个不同的小球放在\(m\)个相同箱子里的方案数（箱子不允许为空）。
-
递推式：\(S(n,m)=S(n-1,m)\times m + S(n-1,m-1)\)，边界为：\(S(0,0)=1\)
证明：
当新加入一个元素时：
- 若将新元素单独放入一个子集，则有\(S(n-1,m-1)\)种方案；
- 若将新元素放入一个现有的非空子集，则有\(S(n-1,m)\)种方案。
根据加法原理，即可得到递推式。
-
一些思考：
- 假如箱子允许为空怎么办？
- 假如是\(m\)个不同的箱子怎么办？
Ans1：考虑剩余箱子个数可能为\([1,m]\)，所以答案为\(\sum_{i=1}^{m}{S(n,i)}\)。
Ans2：考虑到按照第二类斯特林数的方式排列后，\(m\)个箱子的顺序还会对答案造成影响，因此还需要乘上\(m\)个箱子的排列，即\(m!\times S (n,m)\)。
整除分块
-
主要应用场景：在\(O(\sqrt n)\)时间内求出形如\(\sum_{i=1}^{n}{\lfloor \frac{n}{i} \rfloor}\)。
-
大致思路：
整除分块基于这样一个事实：\(\lfloor \frac{n}{i} \rfloor\)的取值总共不会超过\(2\sqrt n\)种。
粗略证明：
- 当\(i \le \sqrt n\)时：由于\(i\)的取值少于\(\sqrt n\)种，所以\(\lfloor \frac{n}{i} \rfloor\)的取值也必然少于\(\sqrt n\)种；
- 当\(i > \sqrt n\)时：\(1 \le \lfloor \frac{n}{i} \rfloor \le \sqrt n\)，所以\(\lfloor \frac{n}{i} \rfloor\)的取值不会超过\(\sqrt n\)种。
故\(\lfloor \frac{n}{i} \rfloor\)不同取值不会超过\(2\sqrt n\)种。
所以，我们可以考虑这样的思路：因为\(\lfloor \frac{n}{i} \rfloor\)的取值是\(\sqrt n\)级别的，所以是需要求出它取每一种值时\(i\)的最大最小值时多少，就可以很快求出答案了。
接下来又需要另外一个结论：如果已知正整数\(p\)，\(n\)（\(n \ge p\)），则使\(\lfloor \frac{n}{i} \rfloor = \lfloor \frac{n}{p} \rfloor\)的最大整数\(i\)的值\(i_{max} = \frac{n}{\lfloor \frac{n}{p} \rfloor}\)。（比较好理解，这里就不证明了）
于是我们只要枚举左端点就好了。
-
模板题：
第二题参考代码：
#include <cstdio> #include <algorithm> #define LL long long using namespace std; const int P = 998244353; LL l,r,sum1,sum2; int main(){ scanf("%lld%lld", &l, &r); l -= 1; for(LL i = 1; i <= l;){ LL j = min(l / (l / i), l); sum1 += (j - i + 1) * (l / i) % P; sum1 %= P; i = j + 1; } for(LL i = 1; i <= r;){ LL j = min(r / (r / i), r); sum2 += (j - i + 1) * (r / i) % P; sum1 %= P; i = j + 1; } printf("%lld\n", (sum2 - sum1 + P) % P); return 0; }
更相损减术
-
可半者半之，不可半者，副置分母、子之数，以少减多，更相减损，求其等也。以等数约之。
翻译：（如果需要对分数进行约分，那么）可以折半的话，就折半（也就是用2来约分）。如果不可以折半的话，那么就比较分母和分子的大小，用大数减去小数，互相减来减去，一直到减数与差相等为止，用这个相等的数字来约分。
-
使用步骤：
第一步：任意给定两个正整数；判断它们是否都是偶数。若是，则用\(2\)约简；若不是则执行第二步。
第二步：以较大的数减较小的数，接着把所得的差与较小的数比较，并以大数减小数。继续这个操作，直到所得的减数和差相等为止。
则第一步中约掉的若干个\(2\)的积与第二步中等数的乘积就是所求的最大公约数。
其中所说的“等数”，就是公约数。求“等数”的办法是“更相减损”法。
注意：辗转相减的复杂度为\(O(n)\)，所以不常用于计算$$gcd$$。
-
重要推论：\(gcd(a,b)=gcd(a,a-b)(a>b)\)。
调和级数
-
\(\sum_{i=1}^{n}{\frac{1}{i}=\ln n+\gamma+\epsilon_n}\)，其中\(\gamma\)为欧拉常熟，约为\(0.5772156649\)，而\(\epsilon_n\)约等于\(\frac{1}{2n}\)，并且随着\(n\)的趋近于正无穷而趋近于\(0\)。
-
应用：
- 某些题目中用于计算时间复杂度（如CSP2020 T2）；
- 计算刚才的式子（废话）。
-
模板题：Luogu6028 算数
-
参考代码：
#include <cstdio> #include <cmath> #define db long double #define LL long long #define beta 0.577216 using namespace std; const int maxn = 1e7 + 10; const int N = 10000000; db s[maxn],Ans; LL n; db S(LL loc){ if(loc <= N) return s[loc]; else return beta + log(loc) + 1.0 / (2 * loc); } db getsum(LL l, LL r){return S(r) - S(l - 1);} int main(){ scanf("%lld", &n); for(int i = 1; i <= N; ++ i) s[i] = s[i - 1] + 1.0 / i; for(LL i = 1; i <= n;){ LL j = n / (n / i); LL val = n / i; Ans += (db)val * getsum(i, j); i = j + 1; } printf("%.10Lf\n", Ans); return 0; }
欧拉函数
-
欧拉函数\(\varphi(n)\)表示\([1,n]\)中与\(n\)互质的数的个数。
-
一些性质：
- 欧拉函数是积性函数。即如果\(gcd(a,b)=1\)，则\(\varphi(a\times b)=\varphi(a)\times \varphi(b)\)。特别地，当\(n\)是奇数时\(\varphi(2n)=\varphi(n)\)。
- \(n=\sum_{d\|n}{\varphi(d)}\)。
- 若\(\varphi(p^k)=p^k-p^{k-1}\)，其中\(p\)为质数。
- 设\(n=\prod_{i=1}^{s}{p_i^{k_i}}\)，其中\(p_i\)为质数，那么\(\varphi(n)=n\times \prod_{i=1}^{s}{\frac{p_i-1}{p_i}}\)。
- \(n（n>1）\)中与\(n\)互质的数的和为\(\frac{\varphi(n)\times n}{2}\)。（由更相损减术推论可以推出）
以上性质的[证明](欧拉函数 - OI Wiki (oi-wiki.org))。
-
如何求欧拉函数值：
-
单个数求值：根据最后一个性质分解质因数就行了。复杂度\(O(\sqrt{n})\)。
-
线性筛：
我们知道，在普通的线性筛中，\(n\)会被最小的质数\(p\)筛去，于是设\(n=n'\times P\)。
分类讨论一下，
-
若\(n'\%P=0\)，那么\(n'\)包含了\(n\)的所有因子，则：
\[\begin{align}\varphi(n)&=n\times \prod_{i=1}^{s}{p_i^{k_i}}\\ &=P\times n' \times \prod_{i=1}^{s}{p_i^{k_i}}\\ &=P\times \varphi(n') \end{align} \] -
若\(n'\%P\neq 0\)，那么\(n'\)与\(P\)互质，则\(\varphi(n)=\varphi(n')\times \varphi(P)\)。
参考代码：
void init(){ memset(phi,0,sizeof(phi)); phi[1] = 1; for(int i = 2; i <= n; ++ i){ if(!vis[i]) p[++cnt] = i, phi[i] = i - 1; for(int j = 1; j <= cnt && (LL)i * p[j] <= n; ++ j){ vis[i * p[j]] = 1; if(i % p[j] == 0){ phi[i * p[j]] = p[j] * phi[i]; break; } else phi[i * p[j]] = phi[p[j]] * phi[i]; } } } -
-
-
主要用于解决有关互质，GCD的问题，如：
欧拉定理
-
若\(\gcd(a,m)=1\)，则\(a^{\varphi(m)} \equiv 1\pmod m\)。（证明主要是自己不会证……)
-
拓展欧拉定理：
\[a^b\equiv \begin{cases}a^{b\:mod\:\varphi(p)},&{gcd(a,p)=1}\\ a^b,&gcd(a,p)\neq1,b<\varphi(p)\\ a^{b\:mod\:\varphi(p)+\varphi(p)}, &gcd(a,p)\neq 1,b\geq\varphi(p) \end{cases}\pmod p \]证明地址同上。
第一类斯特林数
-
也叫斯特林轮换数，\(S(n,k)\)表示将\(n\)个两两不同的元素，划分成\(k\)个非空园排列的方案数。
-
递推式：\(S(n,k)=S(n-1,k-1)+(n-1)\times S(n-1,k)\)。
证明：
与第一类斯特林数类似，考虑新加入一个元素：
- 将该元素放入一个单独的园排列中，共有\(S(n-1,k-1)\)种方案；
- 将该元素插入任意一个已有的园排列中，共有\((n-1)\times S(n-1,k)\)种方案（因为有\(n-1\)个位置可选）。
中国剩余定理（CRT）
-
用于求解线性同余方程组，如：
\[\begin{cases} x &\equiv a_1 \pmod {n_1} \\ x &\equiv a_2 \pmod {n_2} \\ &\vdots \\ x &\equiv a_k \pmod {n_k} \\ \end{cases} \] -
一般步骤：
-
计算所有模数的积\(n\)；
-
对于第\(i\)个方程：
a. 计算\(m_i=\frac{n}{n_i}\)；
b. 计算\(m_i\)在模\(n_i\)意义下的逆元\(m^{-1}\)；
c. 计算\(c_i=m_i\times m_i^{-1}\)（不要对\(n_i\)取模！）。
-
方程组的唯一解为：\(a=\sum_{i=1}^{k}{a_i\times c_i\pmod n}\)。
参考代码：
N = 1; for(int i = 1; i <= n; ++ i) N *= n[i]; for(int i = 1; i <= n; ++ i){ int m = N / n[i]; int _m,y; ex_gcd(m, n[i], _m, y); if(_m < 0) _m += n[i]; Ans += a[i] * _m * m % N; Ans %= N; } -
-
证明：暂无（逃
-
应用：
-
在某些题中，处于某些原因，给出的模数不是质数,但是对其质因数分解以后发现它没有平方因子，那么我们可以对这些质数的质因子进行计算，最后用\(CRT\)整合答案。
例如：SDOI2010古代猪文
-
当\(g=999911659\)时，答案为\(0\)。
-
否则，根据欧拉定理，
\[g^{\sum_{d\|n}{\C^{d}_{n}}}\mod999911659=g^{\sum_{d\|n}{\C^{d}_{n}} \mod 999911658}\mod 999911659 \]所以就是要想办法求出\(\sum_{d\|n}{\C^{d}_{n}}\mod 999911658\)。但是\(999911658\)不是质数，无法直接计算，但是注意到\(999911658=2\times 3\times 4679\times 35617\)，没有平方因子，所以我们可以分别求出\(\sum_{d\|n}{\C^{d}_{n}}\)在模\(2,3,4679,35617\)时的值(使用卢卡斯定理），再用中国剩余定理整合答案。
也就是说，要求一个这样的方程组的解：
\[\begin{cases}x &\equiv a_1 \pmod {2} \\x &\equiv a_2 \pmod {3} \\ x &\equiv a_3 \pmod {4679} \\x &\equiv a_4 \pmod {35617} \\\end{cases} \]参考代码：
#include <cstdio> #include <cstring> #define int long long using namespace std; const int Max = 35617; int Num[Max + 10]; int n,a[20],b[20],N,Ans,num,g; int P; void init(){ memset(Num,0,sizeof(Num)); Num[1] = Num[0] = 1; for(int i = 2; i <= P; ++ i) Num[i] = Num[i - 1] * i % P; } int Pow(int x, int y){ int ans = 1; while(y){ if(y & 1) ans *= x, ans %= P; x *= x, x %= P; y >>= 1; } return ans; } void ex_gcd(int a, int b, int &x, int &y){ if(!b){ x = 1, y = 0; return; } ex_gcd(b, a % b, x, y); int z = x; x = y; y = z - (a / b) * y; } int C(int n, int m){ if(n < m) return 0; return Num[n] * Pow(Num[n - m], P - 2) % P * Pow(Num[m], P - 2) % P; } int Lucas(int n, int m){ if(!m) return 1; return (C(n % P, m % P) * Lucas(n / P, m / P)) % P; } int solve(int n){ Ans = 0; int i; for(i = 1; i * i < n; ++ i) if(n % i == 0){ Ans += (Lucas(n, i) + Lucas(n, n / i)) % P; Ans %= P; } if(i * i == n) Ans += Lucas(n, i), Ans %= P; return Ans; } int CRT(){ N = 1; Ans = 0; for(int i = 1; i <= n; ++ i) N *= a[i]; for(int i = 1; i <= n; ++ i){ int m = N / a[i]; int _m,y; ex_gcd(m, a[i], _m, y); if(_m < 0) _m += a[i]; Ans += b[i] * _m * m % N; Ans %= N; } return Ans; } signed main(){ scanf("%lld%lld", &num, &g); if(g == 999911659) return printf("0\n"), 0; n = 4; a[1] = 2, a[2] = 3, a[3] = 4679, a[4] = 35617; for(int i = 1; i <= 4; ++ i){ P = a[i]; init(); b[i] = solve(num); } P = 999911659; printf("%lld\n", Pow(g, CRT())); return 0; }
-
-
BSGS(大步小步法)
-
用于\(O(\sqrt P)\)时间内求解\(a^x \equiv b \pmod P\)。
其中\(a\)与\(P\)互质。方程的解满足\(0\le x< P\)。注意：只需要互质，不要求\(P\)为质数。
-
算法描述：
令\(x=A\lceil \sqrt P \rceil - B\)，其中\(0\le A,B \le \lceil \sqrt P \rceil\)，则有\(a^{A\left \lceil \sqrt P \right \rceil -B} \equiv b \pmod P\)，进而得到\(a^{A\left \lceil \sqrt P \right \rceil} \equiv ba^B \pmod P\)。
我们已知\(a,b\)，所以可以先算出右边的\(ba^B\)的所有取值，枚举\(B\)，用\(hash/map\)存储，然后计算\(a^{A\left \lceil \sqrt P \right \rceil}\)，枚举\(A\)，寻找是否有与之相等的\(ba^B\)。这样就可以得到所有的\(x\)了。
因为\(A,B\)均小于\(\lceil \sqrt P \rceil\)，所以复杂度为\(O(\sqrt P)\)，如果用\(map\)，则需要加一个\(log\)。
-
ex_BSGS（拓展大步小步法）
用于处理\(a\)，\(P\)不互质的情况。
既然他们不互质，那么就让他们变得互质。
设\(GCD=gcd(a,P)\)。若\(GCD\nmid b\)，则原方程无解，否则我们让两边同时除以\(GCD\)，得到
\[\frac{a}{GCD}\cdot a^{x-1}\equiv \frac{b}{GCD}\pmod{\frac{P}{GCD}} \]如果\(a\)与\(\frac{P}{GCD}\)仍不互质，那么将\(\frac{P}{GCD}\)作为\(P\)再次进行上一步。
令\(D=\prod_{i=1}^k{GCD_i}\)，那么：
\[\frac{a}{D}\cdot a^{x-k}\equiv \frac{b}{D}\pmod{\frac{P}{D}} \]其中\(a\)与\(\frac{P}{D}\)互质，所以可以使用\(BSGS\)求出一个\(x\)，最后只需要将\(x+k\)即可。
注意：在实现过程中，有可能出现\(\frac{a}{D}=\frac{P}{D}\)的情况，这时候\(a^{x-k}\equiv 1\pmod{\frac{P}{D}}\)，所以\(x-k=0\)，将\(k\)作为答案返回即可。
代码实现：
//Luogu P4195 【模板】拓展BSGS //c++11下测评 #include <cstdio> #include <map> #include <unordered_map> #include <cmath> #include <algorithm> #define int long long using namespace std; int n,m,P; unordered_map<int,int> M; int Pow(int x, int y){ int ans = 1; while(y){ if(y & 1) ans *= x, ans %= P; x *= x, x %= P; y >>= 1; } return ans; } int ex_gcd(int a, int b, int &x, int &y){ if(!b){x = 1, y = 0; return a;} int ret = ex_gcd(b, a % b, x, y); int z = x; x = y, y = z - (a / b) * y; return ret; } int gcd(int a, int b){return (!b) ? a : gcd(b, a % b);} int BSGS(int a, int b, int Num){ M.clear(); int sq = ceil(sqrt(P)); int ret = Pow(a, sq), num = 1; int x,y; ex_gcd(Num, P, x, y); x = (x % P + P) % P; b *= x, b %= P; M[b] = 0; for(int i = 1; i <= sq; ++ i){ num *= a; num %= P; M[num * b % P] = i; } num = 1; for(int i = 1; i <= sq; ++ i){ num *= ret; num %= P; if(M.find(num) != M.end()){ int X = i * sq - M[num]; X %= P; return (X + P) % P; } } return -1; } int ex_BSGS(int a, int b){ if(!b) return 0; int k = 0, num = 1, GCD = 0; while((GCD = gcd(a, P)) > 1){ if(b % GCD) return -1; k ++, b /= GCD, P /= GCD, num = num * (a / GCD) % P; if(num == b) return k; } num = BSGS(a, b, num); if(num != -1) num += k; return num; } signed main(){ while(scanf("%lld%lld%lld", &n, &P, &m) != EOF && (n || P || m)){ int ret = ex_BSGS(n, m); if(ret == -1) printf("No Solution\n"); else printf("%lld\n", ret); } return 0; }
持续更新……
