暂时没有归类。

凸包

一道模板题

首先凸包的定义:

用不严谨的话来讲，给定二维平面上的点集，凸包就是将最外层的点连接起来构成的凸多边形，它能包含点集中所有的点。

（百度百科给的解释）

因此以下就可以算是一个凸包

求法有两个

求法 1

第一个是把凸包分为上下两部分，对于下半部分不难发现斜率是在递增的，因此考虑先给所有点按照从左到右从下到上排序，然后类似于单调栈维护一个斜率。

然后伪代码大概就长这样

1
2
3
4
5
for(i=1;i<=n;i++)
{
	s[++size]=p[i];
	while(size>=3&&getk(s[size-2],s[size])<getk(s[size-2],s[size-1])) s[size-1]=s[size],size--;
}

getk(node  a,node  b)getk(node~~a,node~~b) 是求AB所在直线斜率，特别注意A,B相等的时候返回一个极大值。

对于上半部分凸包，只需要从右向左枚举即可。

求法2

上面的求法不需要高中数学功底因此可见只能解决较简单的问题对于比较毒瘤的问题，这样进一步写下去就很麻烦，所以考虑求法2

就是Graham扫描法

首先前置芝士

1.极角排序

这里只介绍按照atan2的方法。

atan2(x,y) 即返回tanα=xytan \alpha =\frac{x}{y} 的α\alpha同时α∈[−π,π]\alpha \in [-\pi,\pi]

例如我们以左下角E(优先考虑最下)点对这些点极角排序，显然可见a

ABCDEFG的坐标：

7
-1.56 -1.11
1.68 -1.63
2.12 2.29
-2 1
-2.66 -3.87
2.46 -3.39
-4.7 -1.73

排序的结果是这样的（FBCADG)可以理解为从0度逆时针扫一圈。

atan2 （x,y）当y=0时，tanαtan\alpha不存在 （即tanx=1)返回 0 但是x=0的时候返回-1

2.叉积

如图所示p⃗1∗p⃗2\vec p_1 * \vec p_2即为平行四边形面积

公式 p⃗1∗p⃗2=x1y2−x2y1\vec p_1 * \vec p_2 = x1y2-x2y1

原理可以参考下面这张出，出处大概在这里，由于图片加载太慢我放在图床上面了

同时四边形面积大于0说明p2在p1顺时针方向，如果==0重合，如果小于0就说明在逆时针方向

算叉积的时候可以用这个公式

(bx−ax)∗(cy−ay)−(by−ay)∗(cx−ax)(b_x-a_x)*(c_y-a_y)-(b_y-a_y)*(c_x-a_x)

=(axby−bxay)+(bxcy−cxby)+(cxay−axcy)=(a_xb_y-b_xa_y)+(b_xc_y-c_xb_y)+(c_xa_y-a_xc_y)

=a⃗∗b⃗ +b⃗∗c⃗ +c⃗∗a⃗=\vec a*\vec b~+\vec b*\vec c ~+ \vec c*\vec a

c如果在a右边的话，c⃗∗a⃗\vec c*\vec a 面积就是小于0的

同理，如果上面的式子大于0，就说明c不在凸包里面。

如下图所示

刚开始d,b在栈里面，此时比较(d,c,b)，发现>0则c入栈

同理对于d，比较(b,d,c)也大于0，入栈

运行完毕，发现O→b→c→d→OO\to b \to c \to d \to O就是凸包。

对于小于O的情况的解释：

例如 b,c,N三个点，按照极角排序，就是b,c,N,现在我们假设b,c在栈里面，就要比较(b,N,c),此时小于0，即c并不在凸包上（N取代了他的位置)

实际上推完了上面的东西，你就能写出来代码了.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<queue>
#include<stack>
#include<vector>
#include<map>
#include<set>
#include<algorithm>

using namespace std;

const double eps=1e-3;

int n,m,size;
double ans;

struct node{
	double x,y;
}p[200200],s[200200];

inline int cmp(node a,node b)
{
	double k1,k2;
	k1=atan2((a.y-p[1].y),(a.x-p[1].x));   
	k2=atan2((b.y-p[1].y),(b.x-p[1].x));
	if(k1!=k2) return k1<k2;
	return a.x<b.x;
}

inline double getdis(node a,node b)
{
	return sqrt(1.0*(a.x-b.x)*(a.x-b.x)+1.0*(a.y-b.y)*(a.y-b.y));  //上面讲的很清楚了 
}

inline double getcj(node a,node b,node c)
{
 	return 1ll*(b.x-a.x)*(c.y-a.y)-1ll*(b.y-a.y)*(c.x-a.x);
}
 
inline int bj(node a,node b)
{
	return (a.y>b.y||(a.y==b.y)&&(a.x>b.x));
}

void gettb()
{
 	p[0].x=p[0].y=0x3f3f3f3f;
 	register int i,j;
 	for(i=1;i<=n;i++)
 	{
 		if(bj(p[0],p[i])) p[0]=p[i],j=i;  //找到最小边的点，并且把他放到第一位 
	}
	swap(p[j],p[1]);
	sort(p+2,p+n+1,cmp);  //对后面的点极点排序 
	s[++size]=p[1];s[++size]=p[2];  //让前两个点入栈（他们一定是在凸包上） 
	for(i=3;i<=n;)
	{
		if(size>=2&&getcj(s[size-1],p[i],s[size])>=0) size--;  //就是上面说的 
		else s[++size]=p[i++];
	}
	s[0]=s[size];  //记得凸包要首尾相连 
	for(i=1;i<=size;i++) ans+=getdis(s[i-1],s[i]); //算距离即可 
}

int main()
{
	register int i,j;
	scanf("%d",&n);
	for(i=1;i<=n;i++) scanf("%lf %lf",&p[i].x,&p[i].y);
	gettb();
	printf("%.2lf",ans);
	return 0;
}

旋转卡壳

旋(xuán)转(zhuàn)卡(qia)壳(qiào)

求凸包的直径一种方法

又是一道模板

首先我们求出来一个凸包（绿色），然后才能找直径（红色）（bushi）

其实只要找到最长距离的两个点就可以了

找直径的方式：找两个点，做一条平行线，然后比较距离~~(其实并没有什么用，但是我们要引入旋转卡壳)~~

实际上我们可以转化为两个相邻的点的线段所在的直线到最远的点的距离,这样单纯的算距离就变成了求三角形面积。

我们在做无用功？

不是，于是精彩的来了，（如下图）CB为底的三角形面积是单峰函数，我们扫描出来了CB为底的最大值，比较GC,GB的距离，然后逆时针上去(E→F→B→...E \to F \to B \to ...以此类推)

这时候对底边面积最大的点一定在逆时针移动（用手推一推就好了）

这样不难发现三角形的顶点（这里姑且认为是临边不与这条平行线重合的点）不会再回到原来的点。

因此时间复杂度达到了美妙的O(Nlog2N)O(Nlog_2N)（排序占了大头）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<queue>
#include<stack>
#include<vector>
#include<set>
#include<map>
#include<algorithm>

#define int long long

using namespace std;

struct node{
	int x,y;
}p[200200],s[200200];

int n,size;

inline int cmp(node a,node b)
{
	double k1,k2;
	k1=atan2((a.y-p[1].y),(a.x-p[1].x));
	k2=atan2((b.y-p[1].y),(b.x-p[1].x));
	if(k1==k2) return a.x<b.x;
	return k1<k2;
}

inline int getdis(node a,node b)
{
	return 1ll*(a.x-b.x)*(a.x-b.x)+1ll*(a.y-b.y)*(a.y-b.y);
}

inline int getcj(node a,node b,node c)
{
	return 1ll*(b.x-a.x)*(c.y-a.y)-1ll*(b.y-a.y)*(c.x-a.x);
}

inline int bj(node a,node b)
{
	return (a.y>b.y||(a.y==b.y)&&(a.x>b.x));
}

inline void gettb()  //这一部分的代码注释最好在我的凸包讲解那里看（不过我还是复制粘贴了过来） 
{
	p[0].x=p[0].y=0x3f3f3f3f;
	register int i,j;
	for(i=1;i<=n;i++) if(bj(p[0],p[i])) p[0]=p[i],j=i;  //找到最小边的点，并且把他放到第一位 
	swap(p[j],p[1]);
	sort(p+2,p+n+1,cmp); //对后面的点极点排序 
	s[++size]=p[1]; s[++size]=p[2];//让前两个点入栈（他们一定是在凸包上） 
	for(i=3;i<=n;)
	{
		if(size>=2&&getcj(s[size-1],p[i],s[size])>=0) size--; //就是上面说的 (这里的上图在凸包的那里) 
		else s[++size]=p[i++];
	}
	s[0]=s[size]; //记得凸包要首尾相连 
}

inline int getmmx()
{
	int ans=0;
	if(size==2) return getdis(s[0],s[1]);  //只有两个点的时候不能构成三角形，所以要特判 
	int j=2;
	for(int i=0;i<size;i++) 
	{
		while(getcj(s[i],s[i+1],s[j])<getcj(s[i],s[i+1],s[j+1])) j=(j+1)%size;  //实际上这也是在找面积，由于面积只有一个顶点，所以下一个j变小的时候就退出 
		ans=max(ans,max(getdis(s[i],s[j]),getdis(s[i+1],s[j]))); //找到距离最远的一条线段，并且与当前最大值比较 
	}
	return ans;  //最后就是直径了 
}

signed main()
{
	ios::sync_with_stdio(false);
	register int i,j;
	cin>>n;
	for(i=1;i<=n;i++)
	{
		cin>>p[i].x>>p[i].y;
	}
	gettb();
	cout<<getmmx()<<endl;
	return 0;
}

素数（质数）

O(N∗N)O(N*N)

就是判断xx能不能被[2,x−1][2,x-1]整除

O(N∗N)O(N*\sqrt N)

由于约数成对出现，如果假如x=p1p2,那么p1≤(x),p2≥(x)p1\leq\sqrt(x),p2\geq\sqrt(x)

证明不难。

O(N∗log2log2N)O(N*log_2log_2N)

很优秀，小学生应该都能想出来的算法。

1.找到1[]

2.找到[1,N]中3的倍数。

3.4被标记，跳过。

4.找到[1,N]中5的倍数。

……

1
2
3
4
5
6
for(i=2;i<=n;i++)
{
	if(!f[i])
	for(j=i+i;j<=n;j+=i)
	f[j]=1;
}

主要代码就这么点。

O(N)O(N)

就是一个欧拉筛。

首先欧筛的大意就是用pri数组存素数，用v数组标记（也有用来存最小的

不难发现pri素组（存素数的） 保持单调上升（即2,3,5,…)。

1
2
3
4
5
6
7
8
9
10
11
12
13
inline void getprime()
{
	register int i,j;
	for(i=2;i<=maxn;i++)
	{
		if(!vis[i]) p[++cnt]=i;
		for(j=1;i*p[j]<=maxn&&j<=cnt;j++)
		{
			vis[i*p[j]]=1;
			if(i%p[j]==0) break; //这里下面来说
		}	
	}
}

imod  p[j]i \mod p[j]为00时，此时根据我们枚举的顺序,令$s=p_j*i $ 则$ p_j$一定是最小的因子，显然是因为i含有p[j]这个因子，那么此时i就不是最大的因子，会造成重复的枚举，此时break即可。

GCD

欧几里得求GCD

即最大公约数

求法有很多种，但是一般用欧几里得

gcd(a,b)=gcd(b,a%b)gcd(a,b)=gcd(b,a\%b)

拓展欧几里得

当我们想要ax+by=k∗gcd(a,b),k∈Z,K≠0ax+by=k*gcd(a,b),k\in Z , K \neq 0的时候，这个时候就需要拓展欧几里得（至于为啥是这么多,看下面的证明）

假设我们现在求出了d=gcd(a,b)d=gcd(a,b)，我们先令k=1

那么

ax+by=gcd(a,b)adx+bdy=1ax+by=gcd(a,b) \\ \frac{a}{d}x+\frac{b}{d}y=1

显然这个时候d∣a,d∣bd|a,d|b,根据除法的定义我们可以继续推

a0x+b0y=1a0=q1b+r1b0=q2r1+r2r1=q3r2+r3a_0x+b_0y=1 \\ a_0=q_1b+r_1 \\ b_0=q_2r_1+r_2\\ r_1=q_3r_2+r_3\\

这时候根据gcd,我们可以求出更多的rr

这个时候我们可以用gcd求出更多的r然后找到r1，最后就可以找到一组特解

(a0,b0)(a_0,b_0)

那么很自然就会想出exgcd

ax+by=gcd(a,b)=gcd(b,a%b)=bx+(a%b)yax+by=gcd(a,b)=gcd(b,a\%b)=bx+(a\%b)y

bx+(a%b)y=bx+(a−⌊ab⌋)y=ay+bx−b⌊ab⌋y=ay+b(x−⌊ab⌋y)bx+(a\%b)y=bx+(a-\lfloor\frac{a}{b}\rfloor)y=ay+bx-b\lfloor\frac{a}{b}\rfloor y\\=ay+b(x-\lfloor\frac{a}{b}\rfloor y)

但是我们首先要求到gcd才能接x，因此代码就要先注意一下

1
2
3
4
5
6
7
8
int exgcd(int a,int b,int &x,int &y)
{
	if(b==0) {x=1;y=0; return a;}
	int t=exgcd(b,a%b,x,y),tmp=x;
	x=y;
	y=tmp-(a/b)*y;
	return t;
}

另一种求法

设

s=p1k1∗p2k2∗p3k3∗...pnknt=p1k1′∗p2k2′∗p3k3′∗...pnkn′s={p_1}^{k_1}*{p_2}^{k_2}*{p_3}^{k_3}*...{p_n}^{k_n}\\ t={p_1}^{k'_1}*{p_2}^{k'_2}*{p_3}^{k'_3}*...{p_n}^{k'_n}

显然gcd为$

t=p1min(k1,k1′)∗p2min(k2,k2′)∗p3min(k3,k3′)∗...pnmin(kn,kn′)t={p_1}^{min(k_1,k'_1)}*{p_2}^{min(k_2,k'_2)}*{p_3}^{min(k_3,k'_3)}*...{p_n}^{min(k_n,k'_n)}

有什么用？当出题人想考你有gcd⁡(a,x)=c,lcm(b,x)=d\gcd(a,x)=c,lcm(b,x)=d的x的个数的时候，这东西就有用了。

题目

P5451 [THUPC2018]密码学第三次小作业

大意，给定c1,c2,e1,e2,Nc_1,c_2,e_1,e_2,N求mm

c1,c2c_1,c_2与NN互质

{c1≡me1 mod Nc2≡me2 mod N\begin{cases}c_1 \equiv m^{e_1} ~mod~ N \\ c_2 \equiv m^{e_2}~mod~N\end{cases}

我们现在合并成一个式子

c1∗c2≡me1+e2 mod Nc_1*c_2\equiv m^{e_1+e_2}~mod~N

好像并不能发现什么，来试试c1,c2c_1,c_2来个几次方

从互质这个地方出发

gcd(c1,N)=gcd(c2,N)=1gcd(c1,N)=gcd(c2,N)=1

一想到互质，就可以来解xe1+ye2=1xe_1+ye_2=1这个不定方程

那么就可以求出m了

但是x,yx,y不一定都大于1

c−x=(c−1)xc^{-x}=({c^{-1})}^x

所以我们只需要找到它的模NN意义下的逆元就行了，由于N不是质数，所以用EXGCD

另外还要注意到N很大，所以要不用龟速乘要不就用高精或者快速模

1
2
3
inline ll mul(ll x,ll y,ll p){
  return ((x*y-(ll)(((long double)x*y+0.5)/p)*p)%p+p)%p;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<queue>
#include<stack>
#include<vector>
#include<set>
#include<map>
#include<algorithm>

#define int long long

using namespace std;

int x,y,n;
int c1,c2,e1,e2;

int exgcd(int a,int b,int &x,int &y)
{
	if(b==0) 
	{
		x=1,y=0; return a;
	}
	int tmp=exgcd(b,a%b,x,y);
	int t=x;
	x=y;
	y=t-a/b*y;
	return tmp;
}

int gsc(int x,int y,int modp)
{
	int ans=0;
	while(y)
	{
		if(y&1) ans=(ans+x)%modp;
		x=(x+x)%modp;
		y>>=1;
	}
	return ans;
} 

int ksm(int x,int y,int modp)
{
	int ans=1;
	while(y)
	{
		if(y&1) ans=gsc(x,ans,modp);
		x=gsc(x,x,modp);
		y>>=1;
	}
	return ans;
} 

int inv(int x,int modp)
{
	int a,b;
	exgcd(x,modp,a,b);
	return (a+modp)%modp; 
}

signed main()
{
	ios::sync_with_stdio(false);
	register int i,j;
	int t;
	cin>>t;
	while(t--)
	{
		cin>>c1>>c2>>e1>>e2>>n;
		exgcd(e1,e2,x,y);
		if(x<0) x=-x,c1=inv(c1,n);
		if(y<0) y=-y,c2=inv(c2,n);
		cout<<gsc(ksm(c1,x,n),ksm(c2,y,n),n)<<endl;
	}
}

逆元

当我们想要在除法意义下面取模的话

比如ab(modp)\frac{a}{b} \pmod p

我们不可能a%pb%p\frac{a\%p}{b\%p}

所以我们引入逆元

aa的逆元就记为a−1a^{-1}

并且a∗a−1≡1(modp)a*a^{-1}\equiv 1\pmod p

为了好看，我们记a−1=xa^{-1}=x

ax≡1(modp)ax\equiv 1 \pmod p

诶，这个用exgcdexgcd不就可以求出来了吗？

1
2
3
4
5
6
7
8
int exgcd(int a,int b,int &x,int &y)
{
	if(b==0) {x=1;y=0; return a;}
	int t=exgcd(b,a%b,x,y),tmp=x;
	x=y;
	y=tmp-(a/b)*y;
	return t;
}

xx就是逆元

P2613 【模板】有理数取余

费马小定理

a(p−1)≡1(modp),pa^{(p-1)}\equiv 1\pmod p,p为质数

转化一下

a∗ap−2≡1(modp)a*a^{p-2}\equiv 1 \pmod p

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<queue>
#include<stack>
#include<vector>
#include<set>
#include<map>
#include<algorithm>

#define int long long

using namespace std;

const int modp=19260817;

int a,b;

inline void qread(int &x)
{
	x=0;
	static char c=getchar();
	while(!isdigit(c)) c=getchar();
	while(isdigit(c)) {x=x*10+c-'0'; if(x>=modp) x%=modp; c=getchar();}
	return ;
}

int ksm(int x,int y,int modp)
{
	int ans=1;
	while(y)
	{
		if(y&1) ans=(ans*x)%modp;
		x=(x*x)%modp;
		y>>=1;
	}
	return ans;
}

signed main()
{
	ios::sync_with_stdio(false);
	qread(a); qread(b);
	if(b==0)
	{
		cout<<"Angry!";
		return 0;
	}
	cout<<(a*ksm(b,modp-2,modp))%modp<<endl;
	return 0;
}

中国剩余定理

{x≡1(mod 3)x≡1(mod 5)x≡2(mod 7)\begin{cases}x\equiv 1(mod~3)\\x\equiv 1(mod~5) \\x\equiv 2(mod~7)\end{cases}

这个时候先得到

{x≡1(mod 3)x≡1(mod 5)x≡1(mod 7)\begin{cases}x \equiv 1(mod~3)\\x\equiv 1(mod~5) \\x\equiv 1(mod~7)\end{cases}

考虑特解

{x≡1(mod 3)x≡0(mod 5)x≡0(mod 7)\begin{cases}x \equiv 1(mod~3)\\x\equiv 0(mod~5) \\x\equiv 0(mod~7)\end{cases}

此时显然可以表示为

令

s=3,t=5∗7s=3,t=5*7

M=3∗5∗7M=3*5*7

x=1+sa,x=0+tbx=1+sa,x=0+tb

1+sa=tb,tb−sa=11+sa=tb,tb-sa=1

显然也可以写成

tb+sa=1tb+sa=1

那就可以通过解这个方程得到s,然后ans+=s∗(M/s)∗1ans+=s*(M/s)*1

然后记得给ans取模，并且上面这种解法对任意互质的p1,p2,p3...pnp_1,p_2,p_3...p_n的话，就是中国剩余定理

扩展中国剩余定理

当p不互质的时候，显然存在一组gcd(a,b)≠1gcd(a,b)\neq1，从而无法使用exgcd求解

但是我们可以一直合并，让他最终变成一个同余方程，或者用对每组方程找到p作为lcm来处理，不过比上面难多了

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<queue>
#include<stack>
#include<vector>
#include<set>
#include<map>
#include<algorithm>

#define int long long

using namespace std;

const int maxn=100100;

int b[maxn],p[maxn];
int n,m;
int x,y;

int exgcd(int a,int b,int &x,int &y)
{
	if(b==0) {x=1;y=0;  return a;}
	int t=exgcd(b,a%b,x,y),tmp=x;
	x=y;
	y=tmp-(a/b)*y;
	return t;
}

inline int mul(int x,int y,int modp)
{
	int ans=0;
	while(y)
	{
		if(y&1) ans=(ans+x)%modp;
		x=(x+x)%modp;
		y>>=1; 
	}
	return ans;
}


inline int excrt()
{
	int i,j;
	int ans=p[1];
	int nowtot=b[1];
	int lcm=1;
	for(i=2;i<=n;i++)
	{
		int tb,tp,tc;
		tb=nowtot; tp=b[i];
		tc=(p[i]-ans%tp+tp)%tp;
		int gd=exgcd(tb,tp,x,y);
		tc/=gd;
		x=mul(x,tc,tp/gd);
		ans+=x*nowtot;
		nowtot*=(tp/gd);
		ans=((ans%nowtot)+nowtot)%nowtot;
	}
	return (ans%nowtot+nowtot)%nowtot;
}

signed main()
{
	ios::sync_with_stdio(false);
	register int i,j;
	cin>>n;
	for(i=1;i<=n;i++) cin>>b[i]>>p[i];
	cout<<excrt()<<endl;
}

矩阵

[124578]\left[ \begin{array}{ccc} 1 & 2 \\ 4 & 5 \\ 7 & 8 \\ \end{array} \right]

这就是个2∗32*3矩阵

[123456] \left[ \begin{array}{ccc} 1 & 2 &3 \\ 4 & 5 &6\\ \end{array} \right]

这个是3*2的

这两东西相乘，得到一个2*2的矩阵

1
2
3
4
5
6
7
8
9
10
11
12
13
struct mat{
	
	int p[3][3];
	
	mat(){ memset(p,0,sizeof(p));}
	
	void fst()
	{
		int i;
		for(i=1;i<=2;i++) p[i][i]=1; 
	}
	
}a,b;

知道了性质，就可以拓展一下，比如斐波那契数列

设有这么个矩阵

[f[n−1]f[n−2]] \left[ \begin{array}{ccc} &f[n-1]&\\&f[n-2]&\\\end{array} \right]

我们想把他变成

[f[n]f[n−1]] \left[ \begin{array}{ccc} &f[n]&\\&f[n-1]&\\\end{array} \right]

根据矩阵乘法的定义，显然

[1110] \left[ \begin{array}{ccc} 1&1\\1&0\\\end{array} \right]

就满足

同时乘法满足交换律，意思就是我们吧这个矩阵k次方先求出来，然后再乘

[f[2]f[1]] \left[ \begin{array}{ccc} f[2]\\f[1]\\\end{array} \right]

（显然得到的答案并不是第n项，这地方你自己想办法处理）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<queue>
#include<stack>
#include<vector>
#include<set>
#include<map>
#include<algorithm>

#define int long long

using namespace std;

const int modp=1e9+7;

int n;

struct mat{
	
	int p[3][3];
	
	mat(){ memset(p,0,sizeof(p));}
	
	void fst()
	{
		int i;
		for(i=1;i<=2;i++) p[i][i]=1; 
	}
	
}a,b;

mat operator*(mat a,mat b)
{
	mat ans;
	int i,j,k;
	for(i=1;i<=2;i++) for(j=1;j<=2;j++) for(k=1;k<=2;k++) ans.p[i][j]=(a.p[i][k]*b.p[k][j]+ans.p[i][j])%modp;
	return ans;
}

mat ksm(int x)
{
	mat ans,tmp;
	ans.fst();
	while(x)
	{
		if(x&1) ans=ans*a;
		a=a*a;
		x>>=1;
	}
	return ans;
}

signed main()
{
	ios::sync_with_stdio(false);
	cin>>n;
	if(n<3) { cout<<1<<endl; return 0;}
	a.p[1][1]=1;
	a.p[1][2]=1;
	a.p[2][1]=1;
	b.p[1][1]=b.p[1][2]=1;
	mat rans;
	rans=b*ksm(n-2);
	cout<<rans.p[1][1];
	return 0;
}

Lucas定理

对质数pp:

Cnm≡Cn%pm%p∗Cn/pm/p(mod p)C_n^m\equiv C_{n\%p}^{m\%p}*C_{n/p}^{m/p}(mod ~p)

但是为了方便递归理解，也有写成

Lucasnm≡Lucasn%pm%p∗Cn/pm/p(mod p)Lucas_n^m\equiv Lucas_{n\%p}^{m\%p}*C_{n/p}^{m/p}(mod ~p)

显然这是个简单到不能再简单的结论，所以我们考虑证明

设CabC_a^b，并且用pp进制表示a,ba,b

a=a0+a1p+a2p2+.....+anpna=a_0+a_1p+a_2p^2+.....+a_np^n

b=b0+b1p+b2p2+.....+bnpnb=b_0+b_1p+b_2p^2+.....+b_np^n

现在我们求一下(1+x)p(1+x)^p

其实不用求，显而易见[1,p−1][1,p-1]的系数都包含了p的倍数，那么

(1+x)p≡1+xp(mod p)(1+x)^p\equiv1+x^p (mod~p)

这东西有啥用？

回头来看命题

Cnm≡Cn%pm%p∗Cn/pm/p(mod p)C_n^m\equiv C_{n\%p}^{m\%p}*C_{n/p}^{m/p}(mod ~p)

换成数学上的东西

Cab≡Ca0b0∗Ca1b1∗...∗Canbn(mod p)C_a^b \equiv C_{a_0}^{b_0}*C_{a_1}^{b_1}*...*C_{a_n}^{b_n}(mod ~p)

另外可以发现（1+x）a（1+x）^a可以被拆分成

(1+x)a=(1+x)a0∗(1+x)a1p∗(1+x)a2p2∗...∗((1+x)anpn)(1+x)a≡(1+x)a0∗(1+xp)a1∗(1+xp2)a2∗...∗(1+xpn)an(mod p)(1+x)^a=(1+x)^{a_0}*(1+x)^{a_{1}p}*(1+x)^{a_{2} p_2}*...*((1+x)^{a_{n}p_n})\\ (1+x)^a\equiv(1+x)^{a_0}*(1+x^p)^{a_1}*(1+x^{p^2})^{a_2}*...*(1+x{p^n})^{a_n}(mod~p)

根据二项式定理，CabC_a^b就是xbx^b项的系数，两边都拿出上面那个项，就和c与b有关，两边都是b，自然约掉了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<queue>
#include<stack>
#include<vector>
#include<set>
#include<map>
#include<algorithm>

#define int long long

using namespace std;

int n,m,p;
int a[200200];

int qpow(int x, int y)
{
	int ans=1;
	x%=p;
	while(y)
	{
		if(y&1) ans=(ans*x)%p;
		x=(x*x)%p; 
		y>>=1;
	}
	return ans;
}

int c(int x,int y)
{
	if(x<y) return 0;
	return ((a[x]*(qpow(a[y],p-2))%p)*qpow(a[x-y],p-2))%p;
}

int getlucas(int x,int y)
{
	if(y==0) return 1;
	return (c(x%p,y%p)*getlucas(x/p,y/p))%p;
}

signed main()
{
	int t;
	register int i,j;
	cin>>t;
	a[0]=1;
	while(t--)
	{
		cin>>n>>m>>p;
		for(i=1;i<=p;i++) a[i]=(a[i-1]*i)%p;
		cout<<getlucas(n+m,n)<<endl;
	} 
}

扩展卢卡斯

先给p分解质因数，然后用一个中国剩余定理的方程组合并即可。

但是m!(n−m)!m!(n-m)!不保证不与pkp^k互质，所以都要先除掉它的因子，并且n!n!这种式子在%p\%p情况下会出现循环节，以及一部分可以用来递归求解，建议看洛咕的题解（公式太多懒得打）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<queue>
#include<stack>
#include<vector>
#include<set>
#include<map>
#include<algorithm>

#define int long long

using namespace std;

int n,m,p;
int x,y;

int exgcd(int a,int b,int &x,int &y)
{
	if(b==0) {x=1; y=0;return a;}
	int tmp=exgcd(b,a%b,x,y),t=x;
	x=y;
	y=t-a/b*y;
	return tmp;
}

inline int getinv(int a,int p)
{
	exgcd(a,p,x,y);
	return ((x%p)+p)%p;
}

int ksm(int x,int y,int modp)
{
	int ans=1;
	while(y)
	{
		if(y&1) ans=(ans*x)%modp;
		x=(x*x)%modp;
		y>>=1;
	}
	return ans;
}

int f(int x,int pi,int pk)
{
	if(x==0) return 1;
	int p1=1;
	int i,j;
	for(i=1;i<=pk;i++) if(i%pi) p1*=i,p1%=pk;
	p1=ksm(p1,x/pk,pk);
	for(i=1;i<=x%pk;i++) if(i%pi) p1*=i,p1%=pk; 
	return (f(x/pi,pi,pk)*p1)%pk;
}

int getk(int x,int p)
{
	int ans=0;
	while(x) ans+=x/p,x/=p;
	return ans; 
}

inline int c(int n,int m,int pi,int pk)
{
	int ta,tb,tc;
	ta=f(n,pi,pk);
	tb=f(m,pi,pk);
	tc=f(n-m,pi,pk);
	int k=getk(n,pi);
	k-=getk(m,pi);
	k-=getk(n-m,pi);
	int ans=((ta*getinv(tb,pk))%pk)*getinv(tc,pk);
	ans=(ans%pk)*ksm(pi,k,pk);
	ans%=pk;
	return ans;
}

inline int crt(int a,int b)
{
	return getinv(p/b,b)%p*(p/b)%p*(a)%p;
}

inline int exlucas(int n,int m)
{
	int i,tmp=p;
	int ans=0,pk=1;
	int len=sqrt(p)+5;
	for(i=2;i<=len;i++)
	{
		if(tmp%i!=0) continue;
		pk=1;
		while(tmp%i==0) pk*=i,tmp/=i;
		ans=(ans+crt(c(n,m,i,pk),pk))%p; 
	}
	if(tmp!=1) ans=(ans+crt(c(n,m,tmp,tmp),tmp))%p; 
	return ((ans%p)+p)%p;
}

signed main()
{
	ios::sync_with_stdio(false);
	cin>>n>>m>>p;
	cout<<exlucas(n,m)<<endl; 
}

1.复数

BSGS

带步小步算法

给定一个质数 pp，以及一个整数 bb，一个整数 nn，现在要求你计算一个最小的 ll，满足bl≡n(modp)b^l \equiv n \pmod p。

设l=k∗a−cl=k*a-c

bk∗a−c≡n(modp)b^{k*a-c}\equiv n \pmod p

bk∗a≡n∗bc(modp)b^{k*a}\equiv n*b^c \pmod p

现在我们需要枚举k和c，所以要找到一个最好的cc使枚举次数最小

用分块思想想一想啊，显然a=pa=\sqrt{p}有最小值。

而且可以保证[0,p−1][0,p-1]为一个循环节。

所以找符合条件的k∗pk*\sqrt{p}就行了

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<queue>
#include<stack>
#include<vector>
#include<set>
#include<map>
#include<algorithm>

#define int long long

using namespace std;

const int max=1001000;

map < int , int > mp;

int p,b,n,m; 

int qpow(int x,int y)
{
	int ans=1;
	while(y)
	{
		if(y&1) ans=(ans*x)%p;
		x=(x*x)%p;
		y>>=1;
	}
	return ans;
}

signed main()
{
	ios::sync_with_stdio(false);
	register int i,j;
	cin>>p>>b>>n;
	int phi=p;
	m=ceil(sqrt(phi));
	j=n;
	for(i=0;i<=m;i++)
	{
		mp[j]=i;
		j=(j*b)%p;
	} 
	int jmp=qpow(b,m);
	j=jmp;
	for(i=1;i<=m;i++)
	{
		if(mp.count(j)) {cout<<i*m-mp[j]<<endl; return 0;}
		j=(j*jmp)%p;
	}
	cout<<"no solution"<<endl;
	return 0;
}

欧拉定理

欧拉函数

ϕ(x)\phi(x)或者φ(x)\varphi(x)表示1→x1\to x中与xx互质的正整数个数

显然ϕ(p)=p−1,p\phi(p)=p-1,p为质数

从这里我们可以推导出

ϕ(pe)=(p−1)∗pe−1\phi(p^e)=(p-1)*p^{e-1}

根据唯一分解定理，有

P=k1e1∗p2e2∗...∗pnenϕ(P)=ϕ(p1e1)∗ϕ(p2e2)∗...∗ϕ(pnen)=Πi=1n(pi−1)∗(piei−1)=Πi=1n(1−1pi)∗(piei)=n∗Πi=1n(1−1pi)\begin{align} P&=k_1^{e_1}*p_2^{e_2}*...*p_n^{e_n}\\ \phi(P)&=\phi(p_1^{e_1})*\phi(p_2^{e_2})*...*\phi(p_n^{e_n})\\ &= \Pi_{i=1}^{n} (p_i-1)*(p_i^{e_i-1}) \\ &= \Pi_{i=1}^{n} (1-\frac{1}{p_i})*(p_i^{e_i}) \\ &= n*\Pi_{i=1}^{n} (1-\frac{1}{p_i}) \end{align}

我们就可以用这种方式求出ϕ(n)\phi(n)

1
2
3
4
5
6
7
8
9
10
void phi_table(int n)
{	
	phi[1]=1;	
	for(int i=2;i<=n;i++) if(!phi[i])//如果i的欧拉函数没有计算则计算	   
	for(int j=i;j<=n;j+=i)//计算素数i及其倍数的欧拉函数，所以进入第二层循环的i都是素因子   
	{	
		if(!phi[j])phi[j]=j;//如果j的欧拉函数没有计算则先赋值为j(由欧拉函数知)
		phi[j]=phi[j]/i*(i-1);//累乘得欧拉函数
	}
}

对于单个n，还有一个O(N)O(\sqrt{N})的方法

1
2
3
4
5
6
7
8
9
10
signed main()
{
	ios::sync_with_stdio(false);
	register int i;
	cin>>m;
	phi=m;
	for(i=2;i*i<=m;i++) 
	if(tmp%i==0){phi-=phi/i; while(tmp%i==0) tmp/=i;}
	if(tmp!=1) phi-=phi/tmp;
}

欧拉定理

aϕ(m)≡1( mod m),a⊥ma^{\phi(m)}\equiv1 (~mod~m),a\perp m

证明：证明这个数有ϕ\phi个不同的同余系，并且一一对应（证明同ap−1a^{p-1}在模pp意义下为aa的逆元）

扩展欧拉定理

先咕一下证明

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<queue>
#include<stack>
#include<vector>
#include<set>
#include<map>
#include<algorithm>

#define int long long

using namespace std;

int a,b,m,phi,tmp,bj;

inline void qread(int &x)
{
	x=0;
	static char c=getchar();
	while(!isdigit(c)) c=getchar();
	while(isdigit(c)) {x=x*10+c-'0'; if(x>=phi) x%=phi,bj=1; c=getchar();}
	return ;
}

inline int ksm(int x,int y,int modp)
{
	int ans=1;
	while(y)
	{
		if(y&1) ans=(ans*x)%modp;
		x=(x*x)%modp;
		y>>=1;
	}
	return ans;
}

signed main()
{
	//ios::sync_with_stdio(false);
	register int i,j;
	cin>>a>>m;
	phi=m; tmp=m;
	for(i=2;i*i<=m;i++) 
	if(tmp%i==0){phi-=phi/i; while(tmp%i==0) tmp/=i;}
	if(tmp!=1) phi-=phi/tmp;
	qread(b);
	cout<<ksm(a,bj?b+phi:b,m);
}