裴蜀定理
拓展欧几里得定理
当已知(a)，(b)时，求一组(x)，(y)使得(a imes x + b imes y = GCD(a, b))
-
解一定存在
-
因为(GCD(a,b) = GCD(b, a \% b))，所以(a imes x + b imes y = GCD(b, a \%b ) = b imes x + (a \% b) imes y = b imes x + (a - lfloor frac ab floor imes b) imes y = a imes y + b imes (x - lfloor frac ab floor imes y))。
所以，我们成功地由(a)与(b)的线性组合转化成了(b)与(a \%b)的线性组合。
说人话，我们可以在求(GCD(a,b))的过程中求出一组(x)和(y)了。
应用：
-
求解线性同余方程
定理1：对于方程(a imes x + b imes y = c)（等价于(a imes x equiv cpmod b)）有整数解的充要条件为，(GCD(a,b)mid c)。（裴蜀定理）
根据定理一，我们可以用拓展欧几里得定理求出(a imes x + b imes y = GCD(a,b))的一组解，再(div GCD(a,b))并( imes c)就行了。
-
求上面那个方程的最小解
定理2：(x_i = x_0 + frac{b}{GCD(a,b)} imes t)。其中，(x_i)为不定方程的任意解，(x_0)为已知解，(t)为任意正整数。
由定理2可知，(x_{min} = x_0 + frac{b}{GCD(a,b)} imes 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}) imes (p_{2} ^ {k_2}) imes (p_3 ^{k_3} * ) imes cdotcdotcdot imes(p_n^{k_n}) ]其中，(p_i)为素数。
-
约数和公式
对于已经分解的整数(A = (p_1 ^ {k_1}) imes (p_{2} ^ {k_2}) imes (p_3 ^{k_3} * ) imes cdotcdotcdot imes(p_n^{k_n}))
(A)的所有因子之和为：
[S = (1+p_1+p_1^2+p_1^3+cdotcdotcdot+p_1^{k_1}) imes(1+p_2+p_2^2+p_2^3+cdotcdotcdot+p_2^{k_2}) imescdotcdotcdot imes(1+p_n+p_n^2+p_n^3+cdotcdotcdot+p_n^{k_n}) ] -
等比数列求和(用递归二分求(1+p+p^2+p^3+cdotcdotcdot+p^n))
-
若(n)为(0)，则答案为(1)。
-
若(n)为奇数，一共有偶数项，则
[1+p+p^1+p^2+p^3+cdotcdotcdot+p^n\=(1+p^{n/2+1})+p imes(1+p^{n/2+1})+cdotcdotcdot+p^{n/2} imes(1+p^{n/2+1})\=(1+p+p^2+cdotcdotcdot+p^{n/2}) imes(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
", Ans);
return 0;
}
做法源于数学一本通（不然我也想不出这么奇怪的命名方式……）
如何(nlog 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(nlog 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) imes 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! imes S (n,m))。
整除分块
-
主要应用场景：在(O(sqrt n))时间内求出形如(sum_{i=1}^{n}{lfloor frac{n}{i} floor})。
-
大致思路：
整除分块基于这样一个事实：(lfloor frac{n}{i} floor)的取值总共不会超过(2sqrt n)种。
粗略证明：
- 当(i le sqrt n)时：由于(i)的取值少于(sqrt n)种，所以(lfloor frac{n}{i} floor)的取值也必然少于(sqrt n)种；
- 当(i > sqrt n)时：(1 le lfloor frac{n}{i} floor le sqrt n)，所以(lfloor frac{n}{i} floor)的取值不会超过(sqrt n)种。
故(lfloor frac{n}{i} floor)不同取值不会超过(2sqrt n)种。
所以，我们可以考虑这样的思路：因为(lfloor frac{n}{i} floor)的取值是(sqrt n)级别的，所以是需要求出它取每一种值时(i)的最大最小值时多少，就可以很快求出答案了。
接下来又需要另外一个结论：如果已知正整数(p)，(n)（(n ge p)），则使(lfloor frac{n}{i} floor = lfloor frac{n}{p} floor)的最大整数(i)的值(i_{max} = frac{n}{lfloor frac{n}{p} floor})。（比较好理解，这里就不证明了）
于是我们只要枚举左端点就好了。
-
模板题：
第二题参考代码：
#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 ", (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 ", Ans); return 0; }
欧拉函数
-
欧拉函数(varphi(n))表示([1,n])中与(n)互质的数的个数。
-
一些性质：
- 欧拉函数是积性函数。即如果(gcd(a,b)=1)，则(varphi(a imes b)=varphi(a) imes 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 imes prod_{i=1}^{s}{frac{p_i-1}{p_i}})。
- (n（n>1）)中与(n)互质的数的和为(frac{varphi(n) imes n}{2})。（由更相损减术推论可以推出）
以上性质的[证明](欧拉函数 - OI Wiki (oi-wiki.org))。
-
如何求欧拉函数值：
-
单个数求值：根据最后一个性质分解质因数就行了。复杂度(O(sqrt{n}))。
-
线性筛：
我们知道，在普通的线性筛中，(n)会被最小的质数(p)筛去，于是设(n=n' imes P)。
分类讨论一下，
-
若(n'\%P=0)，那么(n')包含了(n)的所有因子，则：
[egin{align}varphi(n)&=n imes prod_{i=1}^{s}{p_i^{k_i}}\ &=P imes n' imes prod_{i=1}^{s}{p_i^{k_i}}\ &=P imes varphi(n') end{align} ] -
若(n'\%P eq 0)，那么(n')与(P)互质，则(varphi(n)=varphi(n') imes 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 1pmod m)。（证明主要是自己不会证……)
-
拓展欧拉定理：
[a^bequiv egin{cases}a^{b:mod:varphi(p)},&{gcd(a,p)=1}\ a^b,&gcd(a,p) eq1,b<varphi(p)\ a^{b:mod:varphi(p)+varphi(p)}, &gcd(a,p) eq 1,bgeqvarphi(p) end{cases}pmod p ]证明地址同上。
第一类斯特林数
-
也叫斯特林轮换数，(S(n,k))表示将(n)个两两不同的元素，划分成(k)个非空园排列的方案数。
-
递推式：(S(n,k)=S(n-1,k-1)+(n-1) imes S(n-1,k))。
证明：
与第一类斯特林数类似，考虑新加入一个元素：
- 将该元素放入一个单独的园排列中，共有(S(n-1,k-1))种方案；
- 将该元素插入任意一个已有的园排列中，共有((n-1) imes S(n-1,k))种方案（因为有(n-1)个位置可选）。
中国剩余定理（CRT）
-
用于求解线性同余方程组，如：
[egin{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 imes m_i^{-1})（不要对(n_i)取模！）。
-
方程组的唯一解为：(a=sum_{i=1}^{k}{a_i imes c_ipmod 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 imes 3 imes 4679 imes 35617)，没有平方因子，所以我们可以分别求出(sum_{d|n}{C^{d}_{n}})在模(2,3,4679,35617)时的值(使用卢卡斯定理），再用中国剩余定理整合答案。
也就是说，要求一个这样的方程组的解：
[egin{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 "), 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 ", Pow(g, CRT())); return 0; }
-
-
BSGS(大步小步法)
-
用于(O(sqrt P))时间内求解(a^x equiv b pmod P)。
其中(a)与(P)互质。方程的解满足(0le x< P)。注意：只需要互质，不要求(P)为质数。
-
算法描述：
令(x=Alceil sqrt P ceil - B)，其中(0le A,B le lceil sqrt P ceil)，则有(a^{Aleft lceil sqrt P ight ceil -B} equiv b pmod P)，进而得到(a^{Aleft lceil sqrt P ight ceil} equiv ba^B pmod P)。
我们已知(a,b)，所以可以先算出右边的(ba^B)的所有取值，枚举(B)，用(hash/map)存储，然后计算(a^{Aleft lceil sqrt P ight ceil})，枚举(A)，寻找是否有与之相等的(ba^B)。这样就可以得到所有的(x)了。
因为(A,B)均小于(lceil sqrt P ceil)，所以复杂度为(O(sqrt P))，如果用(map)，则需要加一个(log)。
-
ex_BSGS（拓展大步小步法）
用于处理(a)，(P)不互质的情况。
既然他们不互质，那么就让他们变得互质。
设(GCD=gcd(a,P))。若(GCD mid 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 1pmod{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 "); else printf("%lld ", ret); } return 0; }
持续更新……