模板整理——数论部分

- 数论

- 快速幂（(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
",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
",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
",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}mod 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
");
	for(int i=2;i<=n;i++)
	{
		exgcd(i,p);
		printf("%d
",(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
");
	for(int i=2;i<=n;i++)
	{
		inv[i]=(-1ll*(p/i)*inv[p%i]%p+p)%p;
		printf("%d
",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
");
	for(int i=2;i<=n;i++) printf("%lld
",ans[i]%p);
	return 0;
}

- 矩阵快速幂（(mathcal O(n^3log 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^3log l))。

如果是有权图，考虑在离 (s) 点距离为 (1sim max w_i) 的位置建立虚点，也就是说把 (s) 拆为 (max w_i) 个点（包括他自己）。

这样的复杂度是 (mathcal O(n^3w^3log 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
",x[i]);
    return 0;
}

- 康拓展开（(mathcal O(nlog n))）以及逆康拓展开

康拓展开是一种全排列的编码运算（实质上是按字典序排序），定义为：

[X=sum_{i=1}^n a_i(i-1)! ]

(a_i) 表示第 (i) 个数在没有用过的数中（我们从最高位开始统计），有几个数比他小，于是有：

[a_iin[0,i-1] ]

这个过程用树状数组优化即可 (mathcal O(nlog n)) 的求出一个全排列按字典序的排名。

逆康托展开：

考虑每次对排名数取模 ((i-1)!)，然后可以用一些方法（如线段树上二分）同样可以优化到 (mathcal O(nlog n))。

但阶乘这个东西很大，我们也不知道如何把数据搞成 (mathcal O(nlog 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
",(cal(r)-cal(l-1)+MOD)%MOD);
    return 0;
}

- 中国剩余定理（CRT）（(mathcal O(nlog ext{lcm};a_i))）

- 扩展中国剩余定理（CRT） （(mathcal O(nlog ext{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
",(now.b+now.a)%now.a);
    return 0;
}

如果你要看中国剩余定理的特殊写法，这里有两种：

  1. 直接一股脑得合并（(mathcal O(n^2+nlog 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
",ans);
	return 0;
}
  1. 两两不断合并（(mathcal O(nlog 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
",o.b);
	return 0;
}

- 扩展欧拉定理

欧拉定理
若 (gcd (a,m)=1)，则有：

[a^{varphi(m)}equiv1pmod{m} ]

对于快速幂的优化（扩展欧拉定理）：

[a^bequivegin{cases}a^b&b<varphi(m)\a^{b;mod;varphi(m) +varphi(m)}&bgeqslant varphi(m)end{cases} ]

- Lucas（卢卡斯）定理（(mathcal O(p+Tlog 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
",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
",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
",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+nlog 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
",ans%MOD);
    return 0;
}

拓展：

拉格朗日插值可以在 (mathcal O(k)) 或 (mathcal O(klog k)) 时间内求出:

[sum_{i=1}^n i^k ]

的值。

结论：这是一个 (k+1) 次多项式，我比较屑只会 (mathcal O(klog k))。具体见 CF622F。

- FFT（(mathcal O(nlog n))）

可以优化高精度乘法：

- A*B Problem升级版（FFT快速傅里叶）（(mathcal O(nlog n))）

emmmm，我都忘了，重学还要浪费时间。

所以就咕了。

【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/tlx-blog/p/13899764.html