数论模板
数论模板
快速幂
LL power(LL x,LL k,LL mod)
{
LL res=1; x%=mod;
while(k) {
if(k&1) res=res*x%mod;
x=x*x%mod; k>>=1;
}
return res%mod;
}
龟速乘
用快速加代替乘法。
在形如 \(a=b*c \space mod\space d\)
\(b*c\) 会爆 long long 但 \(b,c,d \in [long\space long]\) 时用到.
有效解决了取模中间结果爆 long long 的问题。
LL muler(LL x,LL k,LL mod)
{
LL res=0;
while(k) {
if(k&1) res=(res+x)%mod;
x=(x+x)%mod; k>>=1;
}
return res;
}
光速乘
用途同 龟速乘，只不过是 \(\mathcal O(1)\) 的，更快。
方法就是利用 C++ 的特性，
- long double 溢出时会存高位。
- long long 溢出时会存低位。
具体见蓝书
LL muler(LL x,LL y,LL MOD)
{
LL high=(long double)x/MOD*y;
LL low=x*y-high*MOD;
return (low%MOD+MOD)%MOD;
}
埃氏筛法
typedef long long LL;
const int N=1e6+5;
LL p[N],tot=0;
bool tag[N];
void prime(int n)
{
LL i,j;
for(i=2;i<=n;i++) {
if(tag[i]) continue;
p[++tot]=i;
for(j=i*i;j<=n;j+=i)
tag[j]=true;
}
return;
}
线性筛
const int N=1e6+5;
int p[N],tot=0;
bool tag[N];
void prime(int n)
{
int i,j;
for(i=2;i<=n;i++) {
if(!tag[i]) p[++tot]=i;
for(j=1;j<=tot && p[j]*i<=n;j++) {
tag[p[j]*i]=true;
if(i%p[j]==0) break;
}
}
return;
}
线性筛求欧拉函数
const int N=1e6+5;
int phi[N],p[N],tot=0;
bool tag[N];
void prime(int n)
{
int i,j;
phi[1]=1; // 1 与 1 互质
for(i=2;i<=n;i++) {
if(!tag[i]) p[++tot]=i,phi[i]=i-1; // i 是素数 [1,i-1] in Z 与 i 互质。
for(j=1;j<=tot && p[j]*i<=n;j++) {
tag[p[j]*i]=true;
if(i%p[j]) phi[i*p[j]]=phi[i]*(p[j]-1);
else {
phi[i*p[j]]=phi[i]*p[j];
break;
}
}
}
return;
}
质因数分解
vector< pair<LL,int> > factor;
void calc(LL x)
{
factor.clear();
for(LL i=2;i*i<=x;i++) {
if(x%i==0) {
factor.push_back(make_pair(i,0));
while(x%i==0) {
x/=i;
factor[factor.size()-1].second++;
}
}
}
if(x>1) factor.push_back(make_pair(x,1));
return;
}
时间复杂度：\(O(\sqrt{n})\)
- AcWing867. 分解质因数
如果把质数筛出来，分解更快。
void calc(int d)
{
v.clear();
for(int i=1;i<=tot&&p[i]*p[i]<=d;i++) {
if(d%p[i]==0) {
v.push_back(make_pair(p[i],0));
while(d%p[i]==0) {
d/=p[i];
v[v.size()-1].second++;
}
}
}
if(d>1) v.push_back(make_pair(d,1));
return;
}
时间复杂度约为：\(\mathcal O(\sqrt{n \over \ln n})\),会快上10倍左右。
如果要分解的 \(n \leq 10^7\), 那可以在 线性筛 过程中存下最小质因子，每次把最小质因子加入集合，然后除去，循环往复。
这样可以做到 \(O(\log n)\)
约数分解
1 .试除法求约数
时间复杂度：\(O(\sqrt{n})\)
const int N=1e5+5;
int p[N],tot;
void divide(int x)
{
tot=0; int i;
for(i=1;i*i<=x;i++) {
if(x%i==0) {
p[++tot]=x/i;
p[++tot]=i;
}
}
if((i-1)*(i-1)==x)
tot--;
sort(p+1,p+tot+1);
return;
}
- AcWing.869 试除法求约数
2 . 质因数+dfs求约数
// 质因数分解
void calc(int d)
{
v.clear();
for(int i=1;i<=tot&&p[i]*p[i]<=d;i++) {
if(d%p[i]==0) {
v.push_back(make_pair(p[i],0));
while(d%p[i]==0) {
d/=p[i];
v[v.size()-1].second++;
}
}
}
if(d>1) v.push_back(make_pair(d,1));
return;
}
//用质因数凑约数
vector<LL> factor;
void dfs(int step,LL gma)
{
LL i,j;
if(step==(int)v.size()) {
factor.push_back(gma);
return;
}
for(i=0,j=1;i<=v[step].second;i++,j*=v[step].first)
dfs(step+1,gma*j);
return;
}
主函数部分
calc(x);
dfs(0,1);
时间复杂度约为：\(O(\sqrt{n}/INn)\),会快上10倍左右。
- AcWing.200 Hankson的趣味题
求单个数的欧拉函数
int varphi(int x)
{
int res=x;
int i;
for(i=2;i*i<=x;i++) {
if(x%i==0) {
while(x%i==0) x=x/i;
res=res/i*(i-1);
}
}
if(x>1) res=res/x*(x-1);
return res;
}
扩展欧几里得|线性同余方程求解
\(ax+by=\gcd(a,b)\)
int exgcd(int a,int b,int &x,int &y)
{
if(b==0) {
x=1; y=0;
return a;
}
int z=exgcd(b,a%b,y,x);
y=y-a/b*x;
return z;
}
假设求出了一组特解 \(x_0,y_0\)，那么
\(x=x_0+k*{b \over \gcd(a,b)},y=y_0-k*{a \over \gcd(a,b)},k \in \Z\) 是解系.
更一般的情况
\(ax+by=c\)
先求 \(ax+by=\gcd(a,b)\)
若 \(\gcd(a,b) \not | \:c\) ，根据贝祖定理无解。
否则 \(ax_0+by_0=\gcd(a_,b)\) 故 \(ax_0{c \over \gcd(a,b)}+by_0{c \over \gcd(a,b)}=c\)
则特解为 \(x_0'=_0{c \over \gcd(a,b)},y_0'=y_0{c \over \gcd(a,b)}\)
通解为 \(x=x_0'+k*{b \over \gcd(a,b)},y=y_0'-k*{a \over \gcd(a,b)},k \in \Z\)
**注意只有 \(x_0,y_0\) 乘了 \({c \over \gcd(a,b)}\),后面没有！ **
中国剩余定理 CRT
wiki
所有模数互质。
typedef long long LL;
const int N=15;
// x = ai (mod bi)
// ami * bmi = 1 (mod bi)
LL A[N],B[N];
LL am[N],bm[N];
void exgcd(LL a,LL b,LL &x,LL &y)
{
if(b==0) {
x=1; y=0;
return;
}
exgcd(b,a%b,y,x);
y=y-a/b*x;
return;
}
LL crt(int n)
{
int i;
LL gm=1,res=0,y;
for(i=1;i<=n;i++) gm*=B[i];
for(i=1;i<=n;i++) {
am[i]=gm/B[i];
exgcd(am[i],B[i],bm[i],y);
res=(res+am[i]*bm[i]%gm*A[i]%gm+gm)%gm;
}
return res;
}
exCRT
处理 CRT 中 \(a_i\) 不互质的情况。
一般使用 exCRT 毕竟其又好想又好写。
exCRT 的核心在于维护一个同余方程，并将同余方程不断合并。
假设现在合并了前 \(i-1\) 个方程，不妨设
现在要 \(ans'\) 同时满足两个方程
两式相减有
这是一个同余方程，假设其解为
\(k=y + jp,j \in \Z\)
回代 \(ans'=ans_0+(y+jp)m=(ans_0+ym)+j \cdot pm\)
故 \(ans_0'=ans_0+ym,m'=m \cdot p,ans'=ans_0'+k \cdot m'\)
发现又回到了原来的形式，合并完成。
这个解系满足 \(ans'=(ans_0+ym)+j \cdot pm=ans_0 \pmod m\)
同时也满足第 \(i\) 个方程 $ans'=ans_0+k \cdot m=b_i+k_1 \cdot a_i $
接下来证明这样不会漏解
引理 : 两个同余方程在 \(lcm(a_1,a_2)\) 内有唯一解。
证明 : 反证法。
假设存在两个不同的解 \(0 \le x,y <lcm(a_1,a_2)\)
不妨设 \(x >y\)
则
换句话说
则有
由于 \(x-y <lcm(a_1,a_2)\)
所以 \(x-y=0 ,x=y\) 不符题意。
证毕。
考虑 \(m'=pm\)
因此不会漏解。
在写代码时，直接取 \(m'=lcm(m,a_i)\) 即可。
P4777 【模板】扩展中国剩余定理（EXCRT）
Code：
typedef long long LL;
LL muler(LL x,LL k,LL MOD)
{
LL res=0; k=(k%MOD+MOD)%MOD; x=(x%MOD+MOD)%MOD;
while(k) {
if(k&1) res=(res+x)%MOD;
x=(x+x)%MOD; k>>=1;
}
return res%MOD;
}
LL exgcd(LL a,LL b,LL& x,LL& y)
{
if(b==0) {
x=1; y=0;
return a;
}
LL z=exgcd(b,a%b,y,x);
y-=a/b*x;
return z;
}
LL excrt(int n,LL b[],LL a[])
{
LL m=a[1],ans=b[1]; // 维护方程的解系为 {ans} = ans0 + k*m
for(int i=2;i<=n;i++) { // 尝试把第 i 个方程加入解系。
LL y,z,d=exgcd(m,a[i],y,z);
if((b[i]-ans)%d!=0) return -1; // 无法合并同余方程 --> 没有合法解
y=muler(y,(b[i]-ans)/d,a[i]/d); // 同余方程的最小正整数解。
ans+=y*m;
m=m/d*a[i]; // m = lcm{a[1]...a[i]} , 不要通过同余方程，这样就很方便，注意先除后乘。
ans=(ans%m+m)%m; // 让 ans 为最小正整数。
}
return ans;
}
费马小定理&欧拉定理
注意欧拉定理只是保证存在 \(a^{\varphi(m)}=1 \space \pmod m，\)其中 \(\gcd(a,m)=1\).
\(\varphi(m)\) 既不是唯一使 \(a^x=1 \space \pmod m\)，成立的数，也不是最小的。
但可以证明 最小的 \(x \mid \varphi(m)\),
扩展欧拉定理
\(if (b <m) \: a^b =a^b \pmod m\)
\(if (b \ge m) \: a^b=a^{(b \bmod \varphi(m))+\varphi(m)}\)
BSGS
求解 \(a^x=b \pmod p\),必须满足 \(a,p\) 互质。
const LL INF=0x3f3f3f3f3f3f3f3f;
LL bsgs(LL a,LL b,LL p)
{
if(1%p==b%p) return 0;
unordered_map<LL,LL> mp;
LL k=sqrt(p)+1,ak=1;
for(LL i=0,j=b%p;i<k;i++,j=j*a%p)
mp[j]=i;
for(LL i=0;i<k;i++) ak=ak*a%p;
for(LL i=1,j=ak;i<=k;i++,j=j*ak%p)
if(mp.count(j)) return i*k-mp[j];
return -INF;
}
exBSGS
求解 \(a^x=b \pmod p\), \(a,p\) 不一定互质。
const LL INF=0x3f3f3f3f3f3f3f3f;
LL exgcd(LL a,LL b,LL& x,LL& y)
{
if(b==0) {
x=1,y=0;
return a;
}
LL z=exgcd(b,a%b,y,x);
y-=(a/b)*x;
return z;
}
LL gcd(LL a,LL b)
{
if(b==0) return a;
else return gcd(b,a%b);
}
LL inv(LL a,LL p) // 求 a 在 模p 意义下的逆元，保证 a,p 互质。
{
LL x,y;
exgcd(a,p,x,y);
x=(x%p+p)%p; // 最小正整数解。
return x;
}
LL bsgs(LL a,LL b,LL p)
{
if(1%p==b%p) return 0;
unordered_map<LL,LL> mp;
LL k=sqrt(p)+1,ak=1;
for(LL i=0,j=b%p;i<k;i++,j=j*a%p)
mp[j]=i;
for(LL i=0;i<k;i++) ak=ak*a%p;
for(LL i=1,j=ak;i<=k;i++,j=j*ak%p)
if(mp.count(j)) return i*k-mp[j];
return -INF;
}
LL exbsgs(LL a,LL b,LL p)
{
b=(b%p+p)%p;
if(1%p==b%p) return 0;
LL d=gcd(a,p);
if(d==1) return bsgs(a,b,p);
else if(b%d) return -INF;
else return exbsgs(a%(p/d),b/d*inv(a/d,p/d)%(p/d),p/d)+1;
}
模运算法则|同余性质
-
如果 \(a≡b \pmod m\) 且有 \(c≡d \pmod m\),那么下面的模运算律成立：
- \(a+c≡b+d\pmod m\)
- \(a-c≡b-d\pmod m\)
- \(a\times c≡b \times d\pmod m\)
-
当且仅当 \(c,m\) 互质时，有 \(a \times c ≡ b \times c \pmod m\) 等价于 \(a ≡ b\pmod m\),具体操作是两边同乘 \(c\) 在 意义下 \(\bmod m\) 的逆元。
否则同余方程的两边不可以同除一个数。 -
$a/c ≡ b/c \pmod m $ 等价于 $a ≡ b \pmod {m \times c} $
-
以下用“\(\%m\)”代表 \(\bmod m\)
- \((a+b)\%m=((a\%m)+(b\%m))\%m\)
- \((a-b)\%m=((a\%m)+(b\%m)+m)\%m\)
- \(a \times b \%m=((a\%m)\times (b\%m))\%m\)
- $a/b%m=a%(b \times m)/b $
- 特别地，若 \(m \in \{Primes\}\) 且 \(b\%m \not =0\),有 \(a/b=a \times b^{p-2}\)
- 若 \(m \in \{Primes\}\) 且 \(b\%m=0\)，此时不存在逆元,特殊情况特殊处理。
实践中可以将 \(\mod m\) 化成 \(+k*m,k \in Z\),放入式子里。
逆元
如果一个线性同余方程 \(ax=1 \: (mod \: p)\) ，则 \(x\) 称为 \(a \: mod \: p\) 的逆元，记作 \(a^{-1}\) 。
引理: 当且仅当 \(a,p\) 互质时，才存在 \(a^{-1}\)
证明:
- 当 \(a,p\) 互质时，根据 欧拉定理，有 \(a^{\varphi(p)}=1 \pmod p\). 故 \(a^{\varphi(p)-1} = a^{-1} \pmod p.\)
- 当 \(a,p\) 不互质时，把式子转化成
\(ax+kp=1,k\in Z\),
设 \(d=\gcd(a,b),d>1\), 则根据裴蜀定理 \(ax+kp\) 的值取遍 \(\{nd\}\).
而 \(1 \not \in \{nd\}\) , 故此方程无解。
费马小定理：当 \(p\) 是质数且 \(a\) 不是 \(p\) 的倍数时，有
也即
所以逆元为 \(a^{p-2}\).
在此条件下可用 费马小定理 求逆元。
LL inv(LL a,LL p)
{
return power(a,p-2,p);
}
当条件只有 \(a,p\) 互质时，可以通过求解次同余方程得到逆元 \(ax+kp=1\)。
通过上述分析，其必然有解 ( \(\gcd(a,p)=1 \: | \: 1\)).
LL inv(LL a,LL p) // 求 a 在 模p 意义下的逆元，保证 a,p 互质。
{
LL x,y;
exgcd(a,p,x,y);
x=(x%p+p)%p; // 最小正整数解。
return x;
}
此法适用范围更广。
组合数计算
容斥原理
参考了214. Devu和鲜花
const int S=(1<<n);
LL ans=0;
for(i=0;i<S;i++) {
sign=1, x=m+n-1;
for(j=0;(i>>j);j++)
if((i>>j)&1)
sign*=-1,x-=(a[j+1]+1);
ans=(ans+sign*C(x,n-1))%MOD;
}
Problems
- AcWing 214. Devu和鲜花
- AcWing 215. 破译密码
- #20055. 「NOIP2020 模拟赛 B 组 Day2」互质数对
线性筛求Mobius函数
Mobius 函数
$ \mu(n)= \begin{cases} 1,\quad n=1\ (-1)^{s}, \quad n=p_{1}p_{2}...p_{s} \0, \quad otherwise \end{cases} $
通俗一点讲，就是一个数
- 有两个以上的同一个质因子时，Mobius 函数为零，
- 否则有奇数个质因子时为-1，
- 有偶数个质因子时为1，
Mobius 函数是为了解决这一类问题而产生的。
例如询问一堆数 \(S\) 中与 \(30\) 互质的数的数目，
（约定 \(S_i\) 为\(S\) 中有 \(i\) 为约数的数的数目）
\(ans=S_1-S_2-S_3-S_5+S_{6}+S_{10}+S_{15}-S_{30}.\)
Mobius 函数 \(Mobius(i)\) 即为 \(S_{i}\) 的符号，若没有出现为0.
const int N=5e4+5;
int p[N],mobius[N],tot=0;
bool tag[N];
void prime(int n)
{
mobius[1]=1;
int i,j;
for(i=2;i<=n;i++) {
if(!tag[i])
p[++tot]=i,mobius[i]=-1;
for(j=1;j<=tot && p[j]*i<=n;j++) {
tag[p[j]*i]=true;
if(i%p[j]==0) {
mobius[p[j]*i]=0;
break;
}
else mobius[p[j]*i]=mobius[i]*(-1);
}
}
return;
}
数论分块
也称整除分块。
AcWing.199 余数之和
$n \times k - \sum_{i=1}^{n}\lfloor k/i \rfloor \times i $
for(i=1;i<=min(n,k);i=j+1ll) {
t=k/i; j=min(k/t,n);
sum+=(i+j)*(j-i+1ll)/2ll*t;
}
cout<<n*k-sum;
扩展版：同时维护两个数的数论分块
AcWing.215 破译密码
\(\sum _{i=1}^{min(a,b)} \lfloor a/i \rfloor * \lfloor b/i \rfloor * mobius[i]\)
n=min(a,b);
for(i=1;i<=n;i=j+1) {
j=min(n,min(a/(a/i),b/(b/i)));
ans+=(LL)(a/i)*(b/i)*(s[j]-s[i-1]);
}
最大公约数与最小公倍数
最大公约数
LL gcd(LL a,LL b)
{
if(b==0) return a;
else return gcd(b,a%b);
}
令 \(\gcd(a,b)=d\). 有如下结论：
- \(\gcd(a,b)=\gcd(b,a)=d\).
- \(\gcd(\gcd(a,b),c)=\gcd(a,\gcd(b,c))=\gcd(\gcd(a,c),b)\).
- 乘法：\(\gcd(a*c,b*c)=d*c\).
- 除法：若 \(k|a,k|b\) ,则 \(\gcd(a/k,b/k)=d/k\).
- \(\gcd(a,b)=\gcd(b,a-b)=d,\) 其中 \(a \geq b\).
- \(\gcd(a,b)=\gcd(a \bmod b,b)=\gcd(b\bmod a,a)\).
裴蜀定理
有 \(ax+by=z\) 的 \(z\) 的取值集合为 \(\{x|x=kd,d=\gcd(a,b)\},a,b,x,y\in Z\).
其中 \(a,b\) 已知,\(x,y\) 不定。
推论一:使 \(ax+by=z\) 成立的最小正整数 \(z=d\).
推论二:上述定理可推广到任意多个
有解当且仅当 \(z=k \times \gcd(a_1,a_2,...,a_n)\)。
高斯消元 Gauss
参考文章：
- AcWing 883. 高斯消元解线性方程组
- AcWing 884. 高斯消元解异或线性方程组
- P4783 【模板】矩阵求逆
- P3389 【模板】高斯消元法
int gauss()
{
int r,c;
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;
for(int i=0;i<n+1;i++)
swap(a[r][i],a[t][i]);
double rate=a[r][c]; // !
if(fabs(rate)<eps) continue; // !
for(int i=c;i<n+1;i++) a[r][i]/=rate;
for(int i=r+1;i<n;i++) {
rate=a[i][c];
for(int j=c;j<n+1;j++)
a[i][j]-=rate*a[r][j];
}
r++;
}
if(r<n) {
for(int i=r;i<n;i++)
if(fabs(a[i][n])>eps) return 2; // No Solution !
return 1; // Infinite group solutions !
}
for(int i=n-1;i>=0;i--)
for(int j=i-1;j>=0;j--)
a[j][n]-=a[i][n]*a[j][i];
return 0;
}
数论里的一些小结论
- \([1,n]\) 里有 \(x\) 因子的数的个数为 \(n/x\) 个。
- 阶乘分解
- \(n!\) 中有几个因子 \(x\)
LL get(LL n,LL x)
{
LL res=0;
while(n>=x) res+=n/x,n/=x;
return res;
}
- 若 \(a,b\) 互质，则 \(a,b\) 所不能表示出的数中最大的一个是 \(ab-a-b\)。
