裴蜀定理
拓展欧几里得定理
当已知(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; }
持续更新……