前排预警
本文几乎没有理论(也就是不能当做课件来看), 只有板子, 纯复习用。
由于我很菜, 并且时间不太够, 于是只学了一部分。
像什么 BSGS、ExBSGS、 Miller-Rabin, 还有那个大数质因数分解我都不会(淦)。
一些基本操作都是现写的伪代码, 都标了伪(就是可能出错); 其他的Luogu有模板题的都是过了模板题的代码, 都标了过。
有可能会写错
前言
基础数论的大部分内容是取模意义下的运算, 猜测最初引入这些科技的目的是为了代替高精度算法qwq。
基本操作
- (O(sqrt n))时间判断(n)是否是质数。
理论支撑 : 若 (n) 为合数, 则其至少有一个小于等于(sqrt n)的质因子。
代码(伪):
bool is_pm(int n) {
int tmp=sqrt(n);
for(int i=2; i<=tmp; ++i)
if(n%i == 0) return false;
return true;
}
- (O(sqrt n))时间将(n)分解质因数。
理论支撑:(n)至多有一个大于(sqrt n)的质因子(不是一种)。
代码(伪):
int tot=0, p[_], c[_]; //p[i]保存n的第i个质因子, c[i]保存p[i]的次数
void divide(int n) {
tot = 0;
for(int i=2; i<=sqrt(n); ++i) if(n%i==0) {
//不用担心i是否为质数, 因为如果i是合数且i整除n
//那么i的质因子一定整除n, 所以枚举到合数时不会出现 n%i == 0
p[++tot] = i, c[tot] = 0;
while(n%i==0) ++c[tot], n/=i;
}
if(n>1) p[++tot] = n, c[tot] = 1;
return;
}
由于上面两个操作都是依赖于枚举质数的, 所以可以预处理一张质数表来加速枚举。
由于([1,n])中的质数个数约为(frac{n}{ln(n)}), 所以加速后上面两个操作的时间复杂度都为(O(frac{sqrt n}{ln(sqrt n)}))(约为qwq)。
- (O(log_2 a)) or (O(log_2 b)) 时间求(a、b)的最大公约数。(反正就是(O(log))级别的啦qwq)
理论支撑就是代码里的式子(不会证)
代码(伪)
int gcd(int a, int b) {
return (!b) ? a : gcd(b, a%b);
}
- (O(n))时间筛出([1,n])的素数表(还可以筛欧拉函数, 好像还可以筛莫比乌斯函数)
理论支撑:将筛一个合数(x)的方式限定为用这个数的最小质因子(p)和(frac{x}p)筛这个数, 意即将每个合数被筛的次数限定为一次, 故时间复杂度为(O(n))。
筛素数代码(过):
int m, prime[maxn + 15], v[maxn + 15];
//prime[i]表示第k大素数, v[i]表示数i的最小质因子
void euler(int n) { //筛出[1,n]所有素数
for(int i=2; i<=n; ++i) {
if(!v[i]) v[prime[++m] = i] = i; //i为质数
for(int j=1; j<=m; ++j) {
if(prime[j] > n/i || prime[j] > v[i]) break;
//不筛大于n的, 不以限定形式之外的形式筛
v[prime[j] * i] = prime[j];
}
}
return;
}
exgcd
- (O(log_2 a)) or (O(log_2 b)) 求解 (ax + by = gcd(a,b))的一组解(x、y)。
就是exgcd啦。
理论支撑:是可以考场上推出来的东西所以不写了
还有就是一些细节什么的做做题就可以回想起来了。
青蛙的约会
荒岛野人
代码(过)
void exgcd(int a, int b, int &x, int &y) {
//exgcd可以顺便求gcd(a,b), 且加上这个功能后代码还是一行, 但是我懒得写。
!b ? x=1,y=0 : (exgcd(b, a%b, y, x), y -= x*(a/b));
}
逆元相关
(a)在膜(p)意义下有逆元的充要条件是(a perp p)。
- (O(log_2 a)) or (O(log_2 p)) 求(a)在膜(p)意义下的逆元,条件是(aperp p)。
就是用exgcd解(ax equiv 1 (mod p))
欧拉定理
- (O(sqrt n))时间求 (varphi(n))
理论支撑:欧拉函数的公式。
代码(伪):
int phi(int n) {
int res = n;
for(int i=2; i<=sqrt(n); ++i) if(n%i==0) {
res = res / i * (i-1);
while(n%i==0) n/=i;
}
if(n>1) res = res / i * (i-1);
return res;
}
- (O(n))求([1,n])的欧拉函数
理论支撑:欧拉函数的公式(自己推推)
代码(过):
int m, prime[maxn], v[maxn], phi[maxn];
//phi[i]为i的欧拉函数
//其他参考上面的筛素数代码
void euler(int n) {
for(int i=2; i<=n; ++i) {
if(!v[i]) {
v[prime[++m] = i] = i;
phi[i] = i-1;
}
for(int j=1; j<=m; ++j) {
if(prime[j] > n/i || prime[j] > v[i]) break;
v[prime[j] * i] = prime[j];
phi[prime[j] * i] = phi[i] * (i%prime[j] ? prime[j]-1 : prime[j]);
}
}
}
- (O(n))求([1,n-1])在膜(n)意义下的逆元, 条件是模数必须是质数。
理论支撑:
设模数为(m), 则对于一个数(i),
由于(r<i), 故通过这种方法可以线性递推([1,m-1])的逆元。
(再次强调此方法应用的条件是模数为质数)
代码(过):
int n,p;//n<p, p为模数
long long inv[maxn];//inv[i]为i在膜p意义下的逆元
void solve() {
inv[1]=1;
for(int i=2;i<=n;++i) inv[i]=(p-p/i)*inv[p%i]%p;
}
求解同余方程组
就是求
的一个解(x)(通常是求最小正整数解)。
模数两两互质的情况。
这种情况下一定有解且求解比较简单。
这种情况下可以为每个(i)构造一个数(node_i), 满足
且
则(x)的一个可行解就是(sum_{i=1}^n a_i *node_i),另外,(x)的通解是
如何构造(node_i)?
设
则
代码(过):
#define li long long
int n;
li a[maxn], m[maxn];
//a、m意义见上文
li CRT() {
li M = 1ll, res = 0ll;
// res就是要求解的x
for(int i=1; i<=n; ++i) M *= m[i]; //这就是M的意义
for(int i=1; i<=n; ++i) {
li Mi = M/m[i];
li Ii = inv(Mi, m[i]);// inv(a, b)是a在膜b意义下的逆元
res += Mi * Ii % M * a[i] % M;
res %= M;
}
return ((res%M+M)%M);
}
模数不一定两两互质的情况
设(X_k)为前(k)个方程组成的方程组的一个特解, 则前(k)个方程组成的方程组的通解为
若前(k+1)个方程组成的方程组有解, 则方程
有解。
由此, 可以一步一步推出整个方程组的解(有解的话)。
代码实现中有不少细节, 主要是关于解同余方程时对正负数的处理。
无解的情况就是exgcd无解的情况。
代码(过)(那道板子不需要判断无解):
// 各个变量基本都和上面的一样
int n;
li a[maxn], m[maxn];
li EXCRT() {
li X = a[1], LCM = m[1];
for(int i=2; i<=n; ++i) {
// k*LCM +m[i]y = c (mod m[i])
li c = ((a[i]-X)%m[i] + m[i])%m[i];
li x, y, d;
exgcd(LCM, m[i], x, y, d);
x = slow_mul(x, c/d, m[i]/d); //x 修正为最小正整数解; slow_mul(a,b,c) = ( (a*b)%c +c ) % c;
X += x*LCM;
LCM = LCM/d * m[i];
X = (X%LCM+LCM)%LCM;
}
return X;
}
计算膜意义下的组合数
就是计算
(p)为质数的情况。
有个定理叫做(Lucas)定理, 其内容为:
代码(过):
//jc[i]是预处理的i的阶乘
li jc[maxp];
li C(li n, li m, li p) {
if(n<m) return 0ll;
return jc[n] * ksm(jc[m],p-2,p) % p * ksm(jc[n-m],p-2,p) % p;
}
li lucas(li n, li m, li p) {
if(!m) return 1ll;
if(n<m) return 0ll;
return lucas(n/p,m/p,p) * C(n%p,m%p,p) % p;
}
(p)不一定为质数的情况
此时可以用(exLucas)求解。
目标: 求
算法过程:
将(a)质因数分解, 即求出(p_1^{c_1}*p_2^{c_2}*cdots*p_k^{c_k} = a)。
分别计算 (ans_j = C_n^m mod p_j^{c_j} (1 leq j leq k))
最后求出
的解(X), (X)即为(C_n^m mod a)。
如何求解(C_n^m mod p^k (p为质数))?
首先, (C_n^m = frac{n!}{(n-m)!m!}), 但((n-m)!m!)中可能有(p)这个因子, 此时((n-m)!m!)在膜(p)意义下没有逆元, 于是就要将((n-m)!m!)中的(p)因子提出来, 最后再乘回去。 但是提出来的(p)因子是(frac{1}{p})形式的, 还是无法计算, 此时就要将(n!)的(p)因子也提出来。
设(n!)含有(p)因子的个数为(num_n), (m!)含有(p)因子的个数为(num_m),((n-m)!)含有(p)因子的个数为(num_{n-m})。
似乎可以证明(p^{num_n-num_m-num_{n-m}} geq 0)。(我不会证也不想证qwq)
于是最后要算的就是
如何求解(frac{n!}{p^{num_n}} mod p^k)?
首先将(n!)展开为
将其中(p)的倍数都提出来
再操作一下
将第一部分扔掉(因为要求的是(frac{n!}{p^{num_n}} mod p^k)),至此可以发现问题被缩小了。
但是如何求([1,n])中与(p)互质的数的乘积(mod p^k)的值呢?
首先(x equiv x + p^k (mod p^k)), 然后就可以(O(p^k))算([1,n])中与(p)互质的数的乘积(mod p^k)的值了。
代码(过)(完整)(原题链接):
#include<bits/stdc++.h>
using namespace std;
#define li long long
li ksm(li a, li b, li p) {
li res = 1ll;
for(;b;b>>=1, a=(a*a)%p)
if(b&1) res = (res*a)%p;
return res%p;
}
void exgcd(li a, li b, li &x, li &y) {
!b ? x=1,y=0 : (exgcd(b,a%b,y,x), y-=x*(a/b));
}
li inv(li a, li p) {
// ax + py = 1;
li x, y;
exgcd(a, p, x, y);
x = ((x%p+p)%p);
return x;
}
li fac(li n, li p, li pk) {
if(!n) return 1ll;
li res = 1ll;
for(li i=0; i<pk; ++i) if(i%p) res = res*i % pk;
res = ksm(res, n/pk, pk);
for(li i=0; i<=n%pk; ++i) if(i%p) res = res*i % pk;
return res * fac(n/p,p,pk) % pk;
}
li C(li n, li m, li p, li pk) {
if(!m) return 0ll;
li f1 = fac(n,p,pk), f2=fac(m,p,pk), f3=fac(n-m,p,pk), cnt = 0ll;
for(li i=n;i;i/=p) cnt += i/p;
for(li i=m;i;i/=p) cnt -= i/p;
for(li i=n-m;i;i/=p) cnt -= i/p;
return f1 * inv(f2,pk) % pk * inv(f3,pk) % pk * ksm(p,cnt,pk) % pk;
}
int tot;
li ri[10001], riri[10001];
li CRT() {
li M = 1ll, res = 0ll;
for(int i=1; i<=tot; ++i) M *= riri[i];
for(int i=1; i<=tot; ++i) {
int nd = M/riri[i];
res += ri[i] * nd % M * inv(nd,riri[i]) % M;
res %= M;
}
return ((res%M+M)%M);
}
li exlucas(li n, li m, li p) {
for(li i=2; i<=sqrt(p); ++i) {
li md = 1ll;
while(p%i==0) p/=i, md*=i;
if(md>1) ri[++tot] = C(n,m,i,md), riri[tot] = md;
}
if(p>1) ri[++tot] = C(n,m,p,p), riri[tot] = p;
return CRT();
}
int main() {
li n, m, p;
cin >> n >> m >> p;
cout << exlucas(n,m,p);
return 0;
}
原根!
从别人课件扒下来的
定义 (g) 是质数 (p) 的原根, 当且仅当 (1,g,g^2, cdots,g^{p-2}) 互不相同 (当然是膜 (p) 意义下的)。
正整数其实也有原根来着
定理 : 任何质数都有原根
原根的个数不一定是唯一的。
如何判定一个数 (g) 是否是 (p) 的原根?
若 (g) 不是 (p) 的原根, 则存在一正整数 (0 le k < p-1) 使得
(g^k equiv 1 (mod p)),
如此,
再由 (gcd(k, p-1) = xk + y(p-1)),
得出
这样, 就可以枚举 (p-1) 的真因子来判断 (g) 是否为 (p) 的原根。
设存在一个 (p-1) 的真因子 (s) 使得 (g^s equiv 1(mod p)), 则 (g) 不是 (p) 的原根, 反之 (g) 是 (p) 的原根。
其实只要判断形如 "((p-1)/一个质数, 这个质数|(p-1))" 形式的因子就好了。
如何找质数 (p) 的原根?
从 2 开始枚举啊
从 (2) 枚举即可, 一般来说找的挺快,不会解释。
代码?
//找质数 p 的原根 g
//和原根表对了几个,应该阔以保证正确,但是 p 太大会炸, 不知道为啥
#include<bits/stdc++.h>
using namespace std;
#define li long long
li ksm(li a, li b, li p) {
li res = 1ll;
for(;b;b>>=1, a=(a*a)%p)
if(b&1) res=(res*a)%p;
return res%p;
}
const int maxn = 1000005;
li p, g;
li tot, P[maxn], C[maxn];
li factot, fac[maxn];
void getfac(li now, li nownum) {
if(now==tot+1) {
fac[++factot] = nownum;
return;
}
for(int i=0; i<=C[now]; ++i) {
getfac(now+1, nownum);
nownum *= P[now];
}
}
bool jd() {
for(int i=1;i<=factot;++i) if(fac[i]!=p-1 && ksm(g,fac[i],p) == 1ll) {
return false;
}
return true;
}
void gtfac(li x) {
for(li i=2; i<=sqrt(x); ++i) if(x%i==0) {
P[++tot] = i, C[tot] = 0;
while(x%i==0) ++C[tot], x/=i;
}
if(x>1) P[++tot] = x, C[tot]=1;
getfac(1,1ll);
}
int main()
{
cin >> p;
gtfac(p-1);
// for(int i=1;i<=factot;++i) cout<<fac[i]<<'
';
// g=3;
// cout<<jd();
cout<<p<<'
';
if(p==3) {
cout<<2;
return 0;
}
for(g=2ll; g<p-1; g+=1ll) if(jd()){
cout << g;
break;
}
return 0;
}
证明:
每做两次辗转相除, gcd/2, 这样的次数最多只有 log w 次。
小凯的疑惑
ax + by = d
由于 gcd(a,b) = 1, 故这个方程一定有整数解
假设 d 已经固定, 则一组解 (x,y) 可以看成平面直角系上的一个点。
根据解 (x,y), 得出通解 (x+kb, y-ka), (k 是整数), 把这个解集里的点都放在直角坐标系上, 它们在一条斜率为负数且斜率一定的直线上, 直线上相邻两个点的坐标差完全一样;且既然斜率是负数, 那么这条直线一定穿过 2、4 象限, 而且由于 d 是正数, 这条直线一定过第一象限。
其实通过前面的分析, 可以看出直线与 y 的交点的纵坐标 h 无限增大的时候, 一定存在一个 h, 使得直线上的解一定在第一象限或 x 轴正半轴或 y 正半轴上, 那么就猜测 d 增大的时候, 直线与 y 轴的交点的纵坐标增大
还是考虑那个 横坐标+b, 纵坐标-a 就从第二象限跳到第四象限的点,由于 d 要尽量大, 所以 x = -1, 跳之后 x = b-1,无法继续考虑。
由于 d 要尽量大, 所以令