模板整理——数论部分
- 数论
- 快速幂（\(\mathcal O(\log n)\)）
#include"iostream"
#include"cstdio"
#include"cmath"
using namespace std;
#define ll long long
ll b,p,k;
ll quickpow(ll a,ll b,ll p)
{
ll ans=1,base=a%p;
while(b)
{
if(b&1) ans=ans*base%p;
base=base*base%p;
b>>=1;
}
return ans%p;
}
int main()
{
scanf("%lld%lld%lld",&b,&p,&k);
return printf("%lld^%lld mod %lld=%lld\n",b,p,k,quickpow(b,p,k)),0;
}
- 求一个数的欧拉函数（\(\varphi(i)\)）值（\(\mathcal O(\sqrt{n})\)）
//有两种写法
int phi(int n)
{
int ans=n;
for(int i=2;i<=sqrt(n);i++)
{
if(n%i==0)
{
ans=ans/i*(i-1);
while(n%i==0) n/=i;
}
}
if(n>=2) ans=ans/n*(n-1);
return ans;
}
int phi(int n)
{
int ans=1,now;
for(int i=1;i<=sqrt(n);i++)
{
now=1;
if(n%i==0)
{
now=i-1,n/=i;
while(n%i==0) now*=i,n/=i;
}
ans*=now;
}
if(n!=1) ans*=n-1;
return ans;
}
- 线性筛（\(\mathcal O(n)\)）
#include"iostream"
#include"cstdio"
#include"cmath"
using namespace std;
int n,q,s;
bool vis[100000005];
int pri[6000005];
int cnt=0;
void get(int maxn)
{
for(int i=2;i<=maxn;i++)
{
if(!vis[i]) pri[++cnt]=i;
for(int j=1;j<=cnt&&pri[j]*i<=maxn;j++)
{
vis[pri[j]*i]=1;
if(i%pri[j]==0) break;
}
}
return;
}
int main()
{
scanf("%d%d",&n,&q);
get(n);
while(q--) scanf("%d",&s),printf("%d\n",pri[s]);
return 0;
}
- 筛欧拉函数（\(\varphi(i)\)）
#include"iostream"
#include"cstdio"
#include"cmath"
using namespace std;
int n,q,s;
bool vis[10000005];
int pri[1000005],phi[10000005];
int cnt=0;
void get(int maxn)
{
for(int i=2;i<=maxn;i++)
{
if(!vis[i]) pri[++cnt]=i,phi[i]=i-1;
for(int j=1;j<=cnt&&pri[j]*i<=maxn;j++)
{
vis[i*pri[j]]=1;
if(i%pri[j]==0)
{
phi[i*pri[j]]=phi[i]*pri[j];
break;
}
phi[i*pri[j]]=phi[i]*(pri[j]-1);
}
}
return;
}
int main()
{
scanf("%d%d",&n,&q);
get(n);
while(q--) scanf("%d",&s),printf("%d\n",pri[s]);
return 0;
}
- 筛因子个数
void get(int n)
{
for(int i=2;i<=n;i++)
{
if(!vis[i]) pri[++cnt]=i,t[i]=2,e[i]=1;
for(int j=1;j<=cnt&&1ll*pri[j]*i<=(ll)n;j++)
{
vis[i*pri[j]]=1;
if(i%pri[j]==0){t[i*pri[j]]=t[i]/(e[i]+1)*(e[i]+2),e[i*pri[j]]=e[i]+1;break;}
else t[i*pri[j]]=t[i]*2,e[i*pri[j]]=1;
}
}
return;
}
- 筛因子和
void get(int maxn)
{
for(int i=2;i<=maxn;i++)
{
if(!vis[i]) pri[++cnt]=i,t[i]=i+1,e[i]=1;
for(int j=1;j<=cnt&&pri[j]*i<=maxn;j++)
{
vis[pri[j]*i]=1;
if(i%pri[j]==0){t[i*pri[j]]=t[i]*pri[j]+e[i],e[i*pri[j]]=e[i];break;}
t[pri[j]*i]=(1+pri[j])*t[i],e[i*pri[j]]=t[i];
}
}
return;
}
- 乘法逆元1
- 乘法逆元2
可能都略微卡常/jk。
- 费马小定理（单次 \(\mathcal O(\log p)\)）
虽然较慢，但是好写，一般是考场首选。
代码就不给了，具体做法就是用快速幂算出 \(a^{p-2}\bmod p\)。
- 扩展欧几里得算法（exgcd）（单次 \(\mathcal O(\log n)\)）
#include"iostream"
#include"cstdio"
#include"cmath"
using namespace std;
int n,p;
int x,y;
void exgcd(int a,int b)
{
if(a==0){y=1,x=0;return;}
exgcd(b%a,a);
int t=x;
x=y-(b/a)*x;
y=t;
}
int main()
{
scanf("%d%d",&n,&p);
printf("1\n");
for(int i=2;i<=n;i++)
{
exgcd(i,p);
printf("%d\n",(x+p)%p);
}
return 0;
}
- 一种递推实现（筛出所有逆元）（\(\mathcal O(n)\)）
#include"iostream"
#include"cstdio"
#include"cmath"
using namespace std;
#define MAXN 3000005
int n,p;
int inv[MAXN];
int main()
{
scanf("%d%d",&n,&p);
inv[1]=1;
printf("1\n");
for(int i=2;i<=n;i++)
{
inv[i]=(-1ll*(p/i)*inv[p%i]%p+p)%p;
printf("%d\n",inv[i]);
}
return 0;
}
- 利用阶乘的另一种递推实现（好像比较好理解，也可以筛出所有逆元）（\(\mathcal O(n)\)）
#include"iostream"
#include"cstdio"
#include"cmath"
using namespace std;
#define MAXN 3000005
#define ll long long
ll fac[MAXN],ans[MAXN];
int p,n;
ll quickpow(ll a,ll b)
{
ll ans=1,base=a;
while(b)
{
if(b&1) ans=base*ans%p;
base=base*base%p;
b>>=1;
}
return ans%p;
}
int main()
{
scanf("%d%d",&n,&p);
fac[0]=1;
for(int i=1;i<=n;i++) fac[i]=fac[i-1]*(ll)i%p;
ll op=quickpow(fac[n],p-2)%p;
for(int i=n;i>=2;i--)
{
ans[i]=fac[i-1]*op%p;
op=op*i%p;
}
printf("1\n");
for(int i=2;i<=n;i++) printf("%lld\n",ans[i]%p);
return 0;
}
- 矩阵快速幂（\(\mathcal O(n^3\log k)\)）
#include"iostream"
#include"cstdio"
#include"cmath"
using namespace std;
#define MAXN 105
#define ll long long
#define MOD 1000000007
struct Mat
{
ll a[MAXN][MAXN];
}e,o;
int n;
ll k;
Mat mul(Mat x,Mat y,int n)
{
Mat c;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
{
c.a[i][j]=0;//注意清空矩阵
for(int k=1;k<=n;k++)
c.a[i][j]=(c.a[i][j]+x.a[i][k]*y.a[k][j]%MOD)%MOD;
}
return c;
}
Mat quickpow(Mat a,ll b,int n)
{
Mat ans=e,base=a;
while(b)
{
if(b&1) ans=mul(base,ans,n);
b>>=1;
base=mul(base,base,n);
}
return ans;
}
int main()
{
scanf("%d%lld",&n,&k);
for(int i=1;i<=n;i++)
{
e.a[i][i]=1;
for(int j=1;j<=n;j++) scanf("%lld",&o.a[i][j]);
}
o=quickpow(o,k,n);
for(int i=1;i<=n;i++)
{
for(int j=1;j<=n;j++) printf("%lld ",o.a[i][j]%MOD);
puts("");
}
return 0;
}
关于邻接矩阵在图上的应用：
给一张 \(n\) 个点的有向无权图，求从 \(s\) 走到 \(t\) 走恰好 \(l\) 步的方案数。
将邻接矩阵算 \(l\) 次幂，得到的 \(e_{s,t}\) 即为答案。
复杂度为 \(\mathcal O(n^3\log l)\)。
如果是有权图，考虑在离 \(s\) 点距离为 \(1\sim \max w_i\) 的位置建立虚点，也就是说把 \(s\) 拆为 \(\max w_i\) 个点（包括他自己）。
这样的复杂度是 \(\mathcal O(n^3w^3\log l)\)。
对于多次询问的矩阵乘法，考虑预处理出转移矩阵的 \(2^k\) 次幂矩阵，然后对于每次询问只需要离线做即可。
当然这种离线是基于原矩阵只有一行或一列的，这样就不需要 \(\mathcal O(n^3)\) 算乘法了。
更精彩的使用矩阵？看这里->矩阵乘法与斐波那契数列
- 高斯消元（\(\mathcal O(n^3)\)）
//这里应该采用的是高斯消元而非高斯—约旦消元
#include"iostream"
#include"cstdio"
#include"cmath"
using namespace std;
#define eps 10e-8
#define MAXN 105
int n;
double a[MAXN][MAXN];
double x[MAXN],now[MAXN];
int main()
{
scanf("%d",&n);
for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++) scanf("%lf",&a[i][j]);
for(int i=1;i<n;i++)
{
for(int j=i+1;j<=n;j++)
{
for(int k=i;k<=n+1;k++) now[k]=a[i][k]*a[j][i];
for(int k=i;k<=n+1;k++) a[j][k]=a[j][k]*a[i][i];
for(int k=i;k<=n+1;k++) a[j][k]=now[k]-a[j][k];
if(a[j][i+1]==0) return puts("No Solution"),0;
for(int k=n+1;k>i;k--) a[j][k]=a[j][k]/a[j][i+1];
}
}
for(int i=n;i>=1;i--)
{
for(int j=n;j>i;j--)
{
a[i][n+1]=a[i][n+1]-a[i][j]*x[j];
}
x[i]=a[i][n+1]/a[i][i];
}
for(int i=1;i<=n;i++) printf("%.2lf\n",x[i]);
return 0;
}
- 康拓展开（\(\mathcal O(n\log n)\)）以及逆康拓展开
康拓展开是一种全排列的编码运算（实质上是按字典序排序），定义为：
\[X=\sum_{i=1}^n a_i(i-1)!
\]
\(a_i\) 表示第 \(i\) 个数在没有用过的数中（我们从最高位开始统计），有几个数比他小，于是有：
\[a_i\in[0,i-1]
\]
这个过程用树状数组优化即可 \(\mathcal O(n\log n)\) 的求出一个全排列按字典序的排名。
逆康托展开：
考虑每次对排名数取模 \((i-1)!\)，然后可以用一些方法（如线段树上二分）同样可以优化到 \(\mathcal O(n\log n)\)。
但阶乘这个东西很大，我们也不知道如何把数据搞成 \(\mathcal O(n\log n)\) 才能过的东西。
- 凯莱（Cayley）定理
点数为 \(n\) 的完全图的有编号生成树有 \(n^{n-2}\) 种。
- 整除分块（\(\mathcal O(\sqrt{n})\)）
#include"iostream"
#include"cstdio"
#include"cmath"
#include"cstring"
using namespace std;
#define ll long long
#define MOD 998244353
ll l,r;
ll cal(ll n)
{
ll ans=0;
for(ll i=1;i<=n;i++)
{
ans=(ans+(n/(n/i)-i+1)%MOD*(n/i)%MOD)%MOD;
i=n/(n/i);
}
return ans;
}
int main()
{
scanf("%lld%lld",&l,&r);
printf("%lld\n",(cal(r)-cal(l-1)+MOD)%MOD);
return 0;
}
- 中国剩余定理（CRT）（\(\mathcal O(n\log \text{lcm}\;a_i)\)）
- 扩展中国剩余定理（CRT） （\(\mathcal O(n\log \text{lcm}\;a_i)\)）
由于扩展中国剩余定理的代码并不比非扩展的难写多少，并且具有普适性，所以考虑都用扩展中国剩余定理。
然而模板题有毒，所以只有 84 分（指第二个）,可以考虑用 __int128。
#include"iostream"
#include"cstdio"
#include"cmath"
using namespace std;
#define ll long long
int n;
struct qu
{
ll a,b;
}op[100005];
ll x,y;
ll gcd(ll a,ll b){return (a%b==0)?b:gcd(b,a%b);}
void exgcd(ll a,ll b)
{
if(b==0){x=1ll,y=0;return;}
exgcd(b,a%b);
ll t=x;
x=y,y=t-(a/b)*y;
}
ll solve(ll a,ll b,ll c)
{
a=a%c,b=b%c;
ll p=gcd(a,c);
if(b%p!=0) return -1;
a=a/p,b=b/p,c=c/p;
exgcd(a,c);
x=(x*b+c)%c;
return x;
}
qu merge(qu s,qu t)
{
ll rt=(t.b-s.b+t.a)%t.a;
rt=solve(s.a,rt,t.a);
qu ans;
if(rt<0)
{
ans.a=ans.b=-1;
return ans;
}
ans.a=s.a/gcd(s.a,t.a)*t.a;
ans.b=(s.a*rt%ans.a+s.b+ans.a)%ans.a;
return ans;
}
int main()
{
scanf("%d",&n);
for(int i=1;i<=n;i++) scanf("%lld%lld",&op[i].a,&op[i].b);
qu now=op[1];
for(int i=2;i<=n;i++)
{
now=merge(now,op[i]);
if(now.a==now.b&&now.a==-1) return puts("No solution!"),0;
}
printf("%lld\n",(now.b+now.a)%now.a);
return 0;
}
如果你要看中国剩余定理的特殊写法，这里有两种：
- 直接一股脑得合并（\(\mathcal O(n^2+n\log n)\)）
#include<iostream>
#include<cstdio>
using namespace std;
typedef long long ll;
int n;
ll a[15],b[15];
ll sum=1,now=1,ans=0,rt=1;
ll x,y;
inline void exgcd(ll a,ll b)
{
if(b==0)
{
x=1,y=0;
return;
}
exgcd(b,a%b);
ll t=x;
x=y;
y=t-a/b*y;
}
int main()
{
scanf("%d",&n);
for(int i=1;i<=n;i++) scanf("%d%d",&a[i],&b[i]);
for(int i=1;i<=n;i++) rt=rt*a[i];
for(int i=1;i<=n;i++)
{
sum=now=1;
for(int j=1;j<=n;j++) if(i!=j) now=now*a[j],sum=sum*a[j]%a[i];
exgcd(sum,a[i]);
x=(x%a[i]+a[i])%a[i];
ans=(ans+b[i]*x*now)%rt;
}
printf("%lld\n",ans);
return 0;
}
- 两两不断合并（\(\mathcal O(n\log n)\)）
#include<iostream>
#include<cstdio>
using namespace std;
typedef long long ll;
int n;
struct node
{
ll a,b;
}e[15];
ll sum=1,now=1,ans=0,rt=1;
ll x,y;
inline void exgcd(ll a,ll b)
{
if(b==0)
{
x=1,y=0;
return;
}
exgcd(b,a%b);
ll t=x;
x=y;
y=t-a/b*y;
}
ll gcd(ll a,ll b){return b==0?a:gcd(b,a%b);}
node merge(node m,node n)
{
node ans;
ans.a=m.a*n.a;
exgcd(n.a%m.a,m.a);
ll l=(x%m.a+m.a)%m.a;
l=l*m.b%ans.a*n.a%ans.a;
exgcd(m.a%n.a,n.a);
ll r=(x%n.a+n.a)%n.a;
r=r*n.b%ans.a*m.a%ans.a;
ans.b=(l+r)%ans.a;
return ans;
}
int main()
{
scanf("%lld",&n);
for(int i=1;i<=n;i++) scanf("%lld%lld",&e[i].a,&e[i].b);
node o=e[1];
for(int i=2;i<=n;i++) o=merge(o,e[i]);
printf("%lld\n",o.b);
return 0;
}
- 扩展欧拉定理
欧拉定理
若 \(\gcd (a,m)=1\)，则有：
\[a^{\varphi(m)}\equiv1\pmod{m}
\]
对于快速幂的优化（扩展欧拉定理）：
\[a^b\equiv\begin{cases}a^b&b<\varphi(m)\\a^{b\;mod\;\varphi(m) +\varphi(m)}&b\geqslant \varphi(m)\end{cases}
\]
- Lucas（卢卡斯）定理（\(\mathcal O(p+T\log p)\)）
- 迭代实现
#include"iostream"
#include"cstdio"
#include"cmath"
using namespace std;
#define read(x) scanf("%d",&x)
#define ll long long
ll fac[200005];
int t,n,m,p;
void init(int maxn){fac[0]=1ll;for(int i=1;i<=maxn;i++) fac[i]=fac[i-1]*i%p;}
ll quickpow(ll a,ll b)
{
ll ans=1ll,base=a;
while(b)
{
if(b&1) ans=ans*base%p;
b>>=1;
base=base*base%p;
}
return ans%p;
}
ll inv(ll a){return quickpow(a,p-2)%p;}
ll C(int x,int y){return (x>=y)?fac[x]*inv(fac[y]*fac[x-y]%p)%p:0;}
ll lucas(int n,int m)
{
ll ans=1;
while(n>0)
{
ans=ans*C(n%p,m%p)%p;
n/=p,m/=p;
}
return ans;
}
int main()
{
read(t);
while(t--)
{
read(n),read(m),read(p);
init(p+1);
printf("%d\n",lucas(n+m,n)%p);
}
return 0;
}
- 递归实现
#include"iostream"
#include"cstdio"
#include"cmath"
using namespace std;
#define read(x) scanf("%d",&x)
#define ll long long
ll fac[200005];
int t,n,m,p;
void init(int maxn){fac[0]=1ll;for(int i=1;i<=maxn;i++) fac[i]=fac[i-1]*i%p;}
ll quickpow(ll a,ll b)
{
ll ans=1ll,base=a;
while(b)
{
if(b&1) ans=ans*base%p;
b>>=1;
base=base*base%p;
}
return ans%p;
}
ll inv(ll a){return quickpow(a,p-2)%p;}
ll C(int x,int y){return (x>=y)?fac[x]*inv(fac[y]*fac[x-y]%p)%p:0;}
ll lucas(int n,int m)
{
if(n<p) return C(n,m)%p;
else return lucas(n/p,m/p)*C(n%p,m%p)%p;
}
int main()
{
read(t);
while(t--)
{
read(n),read(m),read(p);
init(p+1);
printf("%d\n",lucas(n+m,n)%p);
}
return 0;
}
- 大步小步（BSGS）算法（\(\mathcal O(\sqrt{p}\log p)\)）
#include"algorithm"
#include"iostream"
#include"cstring"
#include"cstdio"
#include"cmath"
using namespace std;
#define ll long long
#define MAXN 200005
#define read(x) scanf("%lld",&x)
ll p,a,b;
int h;
struct num
{
ll val;
int id;
}op[MAXN];
bool cmp(num n,num m){if(n.val==m.val) return n.id>m.id;else return n.val<m.val;}
void get()
{
ll now=b;
for(int i=0;i<h;i++)
{
op[i+1].id=i,op[i+1].val=now;
now=now*a%p;
}
sort(op+1,op+h+1,cmp);
return;
}
int find()
{
ll now,mi=1;
for(int i=1;i<=h;i++) mi=mi*a%p;
now=mi;
for(int i=1;i<=h;i++)
{
int l=1,r=h,mid;
while(l<r)
{
mid=(l+r)>>1;
if(op[mid].val<now) l=mid+1;
else r=mid;
}
if(op[l].val==now) return i*h-op[l].id;
now=now*mi%p;
}
return -1;
}
void work()
{
get();
int rt=find();
if(rt<0) puts("no solution");
else printf("%d\n",rt);
}
int main()
{
read(p),read(a),read(b);
a%=p;
if(b==1) return puts("0"),0;
if(b>=p||a==0) return puts("no solution"),0;
h=ceil(sqrt(p));//注意是上取整
work();
memset(op,0,sizeof(op));
return 0;
}
- 拉格朗日（Lagrange）插值（\(\mathcal O(n^2+n\log n)\)）
#include"iostream"
#include"cstdio"
#include"cmath"
#include"algorithm"
using namespace std;
#define MOD 998244353
#define ll long long
int n;
ll k,x[2005],y[2005];
ll ans=0;
ll quickpow(ll a,ll b)
{
ll ans=1ll,base=a;
while(b)
{
if(b&1) ans=ans*base%MOD;
b>>=1;
base=base*base%MOD;
}
return ans%MOD;
}
ll inv(ll a){return quickpow(a,MOD-2)%MOD;}
int main()
{
scanf("%d%lld",&n,&k);
for(int i=1;i<=n;i++) scanf("%lld%lld",&x[i],&y[i]);
for(int i=1;i<=n;i++)
{
ll now=1ll;
for(int j=1;j<=n;j++) if(i^j) now=now*(x[i]-x[j]+MOD)%MOD;
now=inv(now);
for(int j=1;j<=n;j++) if(i^j) now=now*(k-x[j]+MOD)%MOD;
ans=(ans+y[i]*now%MOD)%MOD;
}
printf("%lld\n",ans%MOD);
return 0;
}
拓展：
拉格朗日插值可以在 \(\mathcal O(k)\) 或 \(\mathcal O(k\log k)\) 时间内求出:
\[\sum_{i=1}^n i^k
\]
的值。
结论：这是一个 \(k+1\) 次多项式，我比较屑只会 \(\mathcal O(k\log k)\)。具体见 CF622F。
- FFT（\(\mathcal O(n\log n)\)）
可以优化高精度乘法：
- A*B Problem升级版（FFT快速傅里叶）（\(\mathcal O(n\log n)\)）
emmmm，我都忘了，重学还要浪费时间。
所以就咕了。
