算法笔记——数学相关

---

个人算法笔记：戳我qwq

高精度

• 核心思想：字符串模拟数组，实现大整数（大于long long的范围）运算。
• 应用范围：

1.$longlong$存不下$QAQ$；
2.模拟运算（高精度本来就是模拟）；

• 读入：

    void init(int a[])
    {
        string s;
        cin>>s;
        a[0]=s.lenth();
        for(int i=1;i<=a[0];++i)
        a[i]=s[a[0]-i]-'0';
    }

• 进位，借位处理

加法进位：

c[i]=a[i]+b[i];
if(c[i]>=10)
{
	c[i]%=10;
	++c[i+1];
}

减法借位：

if(a[i]<b[i])
{
	--a[i+1];
	a[i]+=10;
}
c[i]=a[i]-b[i];

乘法进位：

c[i+j-1]=a[i]*b[j]+x+c[i+j-1];
x=c[i+j-1]/10;
c[i+j-1]%=10;

商和余数：视情况而定

• 代码

高精度加法（模板：A+B problem 高精）

#include<cstdio>
#include<cstring>
using namespace std;

char a1[505],b1[505];
int len1,len2,len3;
int a[505],b[505],c[505];

int main()
{
	scanf("%s%s",a1,b1);
	len1=strlen(a1);
	len2=strlen(b1);
	for(int i=0;i<len1;++i)
	a[len1-i]=a1[i]-'0';
	for(int i=0;i<len2;++i)
	b[len2-i]=b1[i]-'0';
	len3=1;
	int x=0;
	while(len3<=len1||len3<=len2)
	{
		c[len3]=a[len3]+b[len3]+x;
		x=c[len3]/10;
		c[len3]%=10;
		len3++;		
	}
	c[len3]=x;
	if(c[len3]==0)len3--;
	for(int i=len3;i>=1;--i)
	printf("%d",c[i]);
	return 0;
}

高精度减法（模板：高精度减法）

#include<cstdio>
#include<cstring>
#define maxn 10005
using namespace std;

char a1[maxn],b1[maxn],m[maxn];
int len1,len2,len3,flag;
int a[maxn],b[maxn],c[maxn];

int main()
{
	scanf("%s%s",a1,b1);
	if(strlen(a1)<strlen(b1)||(strlen(a1)==strlen(b1)&&strcmp(a1,b1)<0))
	{
		flag=1;
		strcpy(m,a1);
		strcpy(a1,b1);
		strcpy(b1,m);
	}//交换数组
	len1=strlen(a1);
	len2=strlen(b1);
	for(int i=0;i<len1;++i)a[len1-i]=a1[i]-'0';
	for(int i=0;i<len2;++i)b[len2-i]=b1[i]-'0';
	len3=1;
	int x=0;
	while(len3<=len1||len3<=len2)
	{
		if(a[len3]<b[len3])
		{
			a[len3]+=10;
			a[len3+1]--;
		}//借位
		c[len3]=a[len3]-b[len3];
		len3++;
	}
	while(c[len3]==0&&len3>1)len3--;
	if(len3==1&&c[1]==0)
	{
		printf("0");
		return 0;
	}
	if(flag)
	printf("-");
	for(int i=len3;i>=1;--i)
	printf("%d",c[i]);
	return 0;
}

高精度乘法（模板：A*B problem）

#include<cstdio>
#include<cstring>
#define maxn 5005
using namespace std;

char a1[maxn],b1[maxn];
int len1,len2,len3;
int a[maxn],b[maxn],c[maxn];

int main()
{
	scanf("%s%s",a1,b1);
	len1=strlen(a1);
	len2=strlen(b1);
	for(int i=0;i<len1;++i)a[len1-i]=a1[i]-'0';
	for(int i=0;i<len2;++i)b[len2-i]=b1[i]-'0';
	for(int i=1;i<=len1;++i)
	{
		int x=0;
		for(int j=1;j<=len2;++j)
		{
			c[i+j-1]+=a[i]*b[j]+x;
			x=c[i+j-1]/10;
			c[i+j-1]%=10;
		}
		c[i+len2]=x;
	}
	len3=len1+len2;
	while(c[len3]==0&&len3>1)len3--;
	for(int i=len3;i>=1;--i)
	printf("%d",c[i]);
	return 0;
}

高精度除法（高精除以低精）（模板：A/B problem）

#include<cstdio>
#include<cstring>
#define maxn 10005
using namespace std;

int len1,len2;
int b,a[maxn],c[maxn];
char a1[maxn];

int main()
{
	scanf("%s%d",a1,&b);
	len1=strlen(a1);
	for(int i=0;i<len1;++i)a[i+1]=a1[i]-'0';
	int x=0;
	for(int i=1;i<=len1;++i)
	{
		c[i]=(x*10+a[i])/b;
		x=(x*10+a[i])%b;
	}
	len2=1;
	while(c[len2]==0&&len2<len1)len2++;
	for(int i=len2;i<=len1;++i)
	printf("%d",c[i]);
	return 0;
}

高精度除法（高精除以高精）（模板：A/B problem II）

#include<cstdio>
#include<cstring>
#include<iostream>
#define maxn 1005
using namespace std;

int a[maxn],b[maxn],c[maxn];

void read(int a[])
{
	string s;
	cin>>s;
	a[0]=s.length();
	for(int i=1;i<=a[0];++i)
	a[i]=s[a[0]-i]-'0';
}//读入

void print(int a[])
{
	if(a[0]==0)
	{
		printf("0\n");
		return;
	}
	for(int i=a[0];i>0;--i)
	printf("%d",a[i]);
	printf("\n");
	return;
}//输出

int cmp(int a[],int b[])
{
	if(a[0]>b[0])return 1;
	if(a[0]<b[0])return -1;
	for(int i=a[0];i>0;--i)
	{
		if(a[i]>b[i])return 1;
		if(a[i]<b[i])return -1;
	}
	return 0;
}//比较数组a,b的大小

void jian(int a[],int b[])
{
	int flag;
	flag=cmp(a,b);
	if(flag==0)
	{
		a[0]=0;
		return;
	}//相等
	if(flag==1)
	{
		for(int i=1;i<=a[0];++i)
		{
			if(a[i]<b[i])
			{
				a[i+1]--;
				a[i]+=10;
			}//借位
			a[i]-=b[i];
		}
		while(a[0]>0&&a[a[0]]==0)a[0]--;
		return;
	}
}//模拟减法

void numcpy(int p[],int q[],int del)
{
	for(int i=1;i<=p[0];++i)
	q[i+del-1]=p[i];
	q[0]=p[0]+del-1;
}//数组复制

void chugao(int a[],int b[],int c[])
{
	int tmp[maxn];
	c[0]=a[0]-b[0]+1;
	for(int i=c[0];i>=1;--i)
	{
		memset(tmp,0,sizeof(tmp));
		numcpy(b,tmp,i);
		while(cmp(a,tmp)>=0)
		{
			c[i]++;
			jian(a,tmp);
		}	
	}
	while(c[0]>0&&c[c[0]]==0)c[0]--;
	return;
}

int main()
{
	read(a);read(b);
	chugao(a,b,c);
	print(c);
	print(a);//输出余数
	return 0;
}

---

乘法逆元

• 概念

联系小学学过的倒数，$x\ast \frac{1}{x}=1$，在$c++$中如果贸然直接除以$x$，可能会出问题。然而如果对于一个数$x$，若存在另一个数$x^{-1}$，使$x\ast x^{-1}$模上一个数$mod$等于1，那么这个数$x^{-1}$就是$x$的逆元。

• 求法

1.递推法：$inv[i]=-(p/i)\ast inv[pmodi]$(p是模数)
2.拓展欧几里得法：
先预处理一遍阶乘，求$i$得逆元即对$i$的阶乘做一遍拓展欧几里得；

• 应用：

1.逆元求组合数；
2.优雅地除以一个数（化除为乘）；

• 代码
  1.递推法:

#include<bits/stdc++.h>
using namespace std;
#define M 30000000
//map<long long,long long>inv;
int n,p;
long long inv[M+5];
void niyuan(int n,int p)
{
	for(int i=2;i<=n;i++)
	{
		inv[i]=-p/i*inv[p%i];
		while(inv[i]<0)
		inv[i]+=p;
	}
}
int main()
{
	scanf("%d%d",&n,&p);
	inv[1]=1;
	niyuan(n,p);
	for(int i=1;i<=n;i++)
	printf("%I64d\n",inv[i]);
	return 0;
}

2.拓展欧几里得法：
预处理阶乘：

a[0]=1
for(int i=1;i<=n+m;i++)
a[i]=a[i-1]*i%mod;//前缀积

求逆元：

void exgcd(ll a,ll b,ll &x,ll &y)
{
	if(b==0)
	{
		x=1;y=0;
		return; 
	}
	ll xx,yy;
	exgcd(b,a%b,xx,yy);
	x=yy;
	y=xx-(a/b)*yy;
}

ll inv(ll sum)
{
	ll x,y;
	exgcd(sum,mod,x,y);
	x=(x%mod+mod)%mod;
	return x; 
}//求逆元,sum为a[i]

3.快速幂求逆元：
主函数内：

ans=inv(i,mod-2);

求逆元：

int inv(int x,int y)
{
	int ans=1;
	while(y)
	{
		if(y&1)ans=(ans*x)%mod;
		y>>=1;
		x=(x*x)%mod;
	}
	ans%=mod;
	return ans;
}//快速幂求逆元

习题：
$emm$好像逆元没有什么裸题吧，不过倒是有一道模板题：【模板】乘法逆元
还有一道神仙题（滑稽）：【模板】乘法逆元2（其实这道题与逆元没有什么关系）

---

排列组合

• 排列

排列的分类：
1.选排列
最普通的排列，从$m$个数选出$n$个数的方案数。
$A_m^n=(m-n)!$
2.错位排列（错排）
从$m$个数选出$n$个数，数$n$不能在第$n$个位置上的方案数。
$f(n)=(n-1)\ast (f(n-1)+f(n-2))$
3.圆排列
从$n$个数选出$r$个数围成一个圈的方案数。
$A_n^r/r$,若$r=n$，则有(n-1)!种。

• 组合
  $C_m^n=\frac{m!}{n!\ast (m-n)!}$
  模板：逆元求组合数

#include<cstdio>
#define ll long long
using namespace std;
ll t,n,m,mod;
ll a[100005];
void exgcd(ll a,ll b,ll &x,ll &y)
{
	if(b==0)
	{
		x=1;y=0;
		return; 
	}
	ll xx,yy;
	exgcd(b,a%b,xx,yy);
	x=yy;
	y=xx-(a/b)*yy;
}
ll inv(ll sum)
{
	ll x,y;
	exgcd(sum,mod,x,y);
	x=(x%mod+mod)%mod;
	return x; 
}//求逆元
ll C(ll x,ll y,ll mod)
{
	if(x<y)
	return 0;
	return a[x]*inv(a[y])%mod*inv(a[x-y])%mod;
}//求组合数
int main()
{
	a[0]=1;
	scanf("%I64d",&t);
	while(t--)
	{
		scanf("%I64d%I64d%I64d",&n,&m,&mod);
		for(int i=1;i<=n+m;i++)
		a[i]=a[i-1]*i%mod;//前缀积
		printf("%I64d\n",C(n+m-1,m,mod));
	}
	return 0;
}
//方法：求C(n-1,m+n-1)%mod

应用：
1.夹棍法；
2.求杨辉三角；
3.二项式定理；

---

二项式定理

应用：杨辉三角
$(x+y{)}^m={\sum }_{i=0}^{i\le m}{(}_m^i)x^i{y}^{m-i}$

---

质数的判定和应用

• 质数的判定方法

1.暴力枚举法：
对于一个数$x$，从2开始枚举到$\sqrt{x}$(一个数最多有$\sqrt{x}$个约数)，若存在$i$，使得$x$%$i$=0，则$x$为合数，否则为质数。注意特判1和2.
若求1~$n$中的所有质数，那么该算法的复杂度为$O(n^2)$。

2.普通筛法：
假如求1~50的质数，一开始如下：

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

划掉1，从2开始枚举，划掉除2以外2的质数：

	2	3		5		7		9
11		13		15		17		19
21		23		25		27		29
31		33		35		37		39
41		43		45		47		49

接下来，划掉除3以外3的倍数：

	2	3		5		7
11		13				17		19
		23		25				29
31				35		37
41		43				47		49

划掉除5以外5的倍数：

	2	3		5		7
11		13				17		19
		23						29
31						37
41		43				47		49

划掉除7以外7的倍数：

	2	3		5		7
11		13				17		19
		23						29
31						37
41		43				47

到这里我们就筛完了。但细心地你会发现一个问题：一个数可能会被筛多次。举个例子，6在筛2的倍数时会被筛掉，在筛3的倍数时同样会被筛掉，这样就增加了此算法的事件复杂度。有没有办法使一个合数只会被筛一次呢？

3.线性筛：
定义两个数组：$mark[i]$表示$i$是否是质数，$prime[i]$表示从小到大找到的质数。若$mark[i]=0$，则吧$i$加入进质数数组里。接下来，无论$i$是否是质数，从1开始枚举$prime[j]$，将$mark[i\ast prime[j]]$标记为合数，若$i$%$prime[j]$=0，说明$i$是$prime[j]$的倍数，跳出循环，保证每一个合数只会被筛一次。
举个例子，当$i$为4时，质数数组有2和3，枚举质数2，标记42=8为合数。因为4是2的倍数，跳出循环。假如不跳出循环，那么下一个被枚举到的合数就是43=12，而12应该在枚举6*2时筛到，多筛了一遍，增加了时间复杂度。

4.Miller-Rabin素性测试：
玄学算法。
费马小定理：若$p$是质数，则$\forall x^{p-1}\equiv 1(modp)$，但不是所有数只要满足$\forall x^{p-1}\equiv 1(modp)$就是质数，所以光有这一个定理还不够;
二次探测定理：若$p$是质数，且$\forall x^2\equiv 1(modp)$，则$(x+1)(x-1)\equiv 0(modp)$，也就是说$x\equiv \pm 1(modp)$。
Miller-Rabin的流程：
设$p$是大于2的奇数，另$p-1=d\ast 2^r$，d为奇数；
随机选择一个$x$，
首先判断费马小定理是否成立：$x^{p-1}\equiv =1(modp)$；
若符合，则由二项式定理，要么$x^d\equiv 1(modp)$，要么存在一个比$r$小的数$i$满足$x^{d\ast 2^i}\equiv -1(modp)$。
单次复杂度为$O(logn)$，多次随机$x$可保证极高正确率。

• 代码

1.线性筛：

#include<bits/stdc++.h>
using namespace std;
int prime[10000005],mark[10000005],c[100005],tot,m,n;
void su()
{
	int i,j;
	for(i=2;i<=n;i++)
	{
		if(mark[i]==0)//是素数。
		{
			tot++;prime[tot]=i;
		}//增加一个素数。
		for(j=1;j<=tot;j++)
		{
			mark[i*prime[j]]=1;
			if(i%prime[j]==0)
			break;//标记倍数为合数，若遇到倍数，则停止。 
		}
	}
}
int main()
{
	int k;
	scanf("%d%d",&n,&m);
	su();
	mark[1]=1; 
	for(k=1;k<=m;k++)
	scanf("%d",&c[k]);
	for(k=1;k<=m;k++)
	if(mark[c[k]]==0)printf("Yes\n");
	else printf("No\n");
	return 0;
}

2.Miller-Rabin

#include<bits/stdc++.h>
using namespace std;
int prime[10000005],mark[10000005],c[100005],tot,m,n;
void su()
{
	int i,j;
	for(i=2;i<=n;i++)
	{
		if(mark[i]==0)//是素数。
		{
			tot++;prime[tot]=i;
		}//增加一个素数。
		for(j=1;j<=tot;j++)
		{
			mark[i*prime[j]]=1;
			if(i%prime[j]==0)
			break;//标记倍数为合数，若遇到倍数，则停止。 
		}
	}
}
int main()
{
	int k;
	scanf("%d%d",&n,&m);
	su();
	mark[1]=1; 
	for(k=1;k<=m;k++)
	scanf("%d",&c[k]);
	for(k=1;k<=m;k++)
	if(mark[c[k]]==0)printf("Yes\n");
	else printf("No\n");
	return 0;
}#include<bits/stdc++.h>
using namespace std;
int prime[10000005],mark[10000005],c[100005],tot,m,n;
void su()
{
	int i,j;
	for(i=2;i<=n;i++)
	{
		if(mark[i]==0)//是素数。
		{
			tot++;prime[tot]=i;
		}//增加一个素数。
		for(j=1;j<=tot;j++)
		{
			mark[i*prime[j]]=1;
			if(i%prime[j]==0)
			break;//标记倍数为合数，若遇到倍数，则停止。 
		}
	}
}
int main()
{
	int k;
	scanf("%d%d",&n,&m);
	su();
	mark[1]=1; 
	for(k=1;k<=m;k++)
	scanf("%d",&c[k]);
	for(k=1;k<=m;k++)
	if(mark[c[k]]==0)printf("Yes\n");
	else printf("No\n");
	return 0;
}

• 应用

大多数时候作为一种辅助算法，帮助解题（如Pollard-Rho）。
没什么裸题，都是很简单的那种。
但是线性筛素数还可以顺带预处理除一些积性函数，如欧拉函数，莫比乌斯函数等等。
一些题（双倍经验）：
1.樱花
2.樱花

---

约数

• 唯一分解定理
  就是分解质因数。
• 一些性质：

1.将$n$分解质因数：$n=\prod p_i^{r_i}$;
2.约数个数：${\sigma }_0(n)=\prod (1+r_i)$;
3.约数和：${\sigma }_1(n)=\prod {\sum }_{k=0}^{r_i}{p}_i^{}$;
4.约数函数：${\sigma }_t(n)={\sum }_{d|n}dt=\prod {\sum }_{k=0}^{r_i}{p}_i^{kt}$;
以上这些东西，我也看不懂。
例题：
1.[AHOI2007]密码箱（数据过水）;
2.[HAOI2007]反素数;
3.[JLOI2014]聪明的燕姿;

• 最大公因数(gcd)

最大公因数的性质：
1.$gcd(a,b)=gcd(a-b,b)$;
2.$gcd(a,b)=gcd(a\%b,b)$;

最大公因数的求法：
1.更相减损术；
2.欧几里得法。

代码：

#include<cstdio>
using namespace std;

int n,m;

int gcd(int x,int y)
{
	return (y==0)?x:gcd(y,x%y);
}

int main()
{
	scanf("%d%d",&n,&m);
	printf("%d",gcd(n,m));
	return 0;
}

例题：
1.NOIp2009 Hankson的趣味题

---

拓展欧几里得

• 应用
  可求解同余方程$ax\equiv m(modb)$，以及不定方程$ax+by=gcd(a,b)$。

• 方法

1.求解同余方程：
转化为$ax=by+m$，移项$ax-by=m$，另$y=-y$，转化为$ax+by=m$，但注意必须$gcd(a,b)|m$，否则无解。

2.求解不定方程：
对于$ax+by=c$，若满足$gcd(x,y)|c$，方程有解；反之则无解。

• 步骤

斐蜀定理：$ax+by=gcd(a,b)$一定有解。
$gcd(a,b)=gcd(b,a\%b)$。
在欧几里得算法的最后一步时，有$b=0$，$gcd(a,b)=a$，因为$1\ast a+0\ast 0=a$，所以此时$ax+by=gcd(a,b)$有一组解$x=1,y=0$。
当$b!=0$时，递归求解方程$bx+(a\%b)y=gcd(b,a\%b)$，直到$b=0$。

• 代码

#include<bits/stdc++.h>
using namespace std;
int a,b,x,y,s;
int gcd(int a,int b)
{
	int r=a%b;
	if(r==0)
	return r;
	a=b;b=r;
	gcd(a,b);
}
void exgcd(int a,int b,int &x,int &y)
{
	if(b==0)
	{
		x=1;
		y=0;
		return;
	}
	int c=a-(a/b)*b,xx,yy;
	exgcd(b,c,xx,yy);
	x=yy;y=xx-(a/b)*yy;
}
int main()
{
	scanf("%d%d",&a,&b);
	s=gcd(a,b);
	exgcd(a,b,x,y);
	printf("%d %d",x,y);
	return 0;
}

• 应用：
  求解一元二次方程或同于方程，有时也会用来求中国剩余定理和拓展中国剩余定理。

• 例题：

1.NOIP2012 同余方程；
2.青蛙的约会（怎么变蓝题了，我刚做的时候不是绿题吗$QAQ$）

---

大步小步算法（BSGS）

• 应用：求解高次同余方程
  $$a^x\equiv b(modc)$$
  其中$gcd(a,c)=1$。
• 由于我也证不来为什么算法是对的，所以我直接粘贴代码了$qwq$：
  看这篇大佬的博客吧：BSGS及扩展BSGS

#include<cstdio>
#include<cmath>
#include<map>
#define ll long long
using namespace std;

ll a,b,p,ans;
map<ll,ll>mp;
map<ll,bool>vis;

ll quick(ll x,ll y)
{
	ll ret=1;
	while(y)
	{
		if(y&1)ret=(ret*x)%p;
		x=(x*x)%p;
		y>>=1;
	}
	return ret;
}

ll bsgs(ll a,ll b)
{
	mp.clear();vis.clear();
	ll q=sqrt(p);
	b%=p;
	for(ll i=0;i<=q;++i)
	{
		ll ret=b*quick(a,i)%p;
		mp[ret]=i;
		vis[ret]=true;
	}
	a=quick(a,q);
	if(a==0)return b==0?1:-1;
	for(ll i=0;i<=q;++i)
	{
		ll ret=quick(a,i);
		if(vis[ret])
		{
			ll j=mp[ret];
			if(j>=0&&i*q-j>=0)
			return i*q-j;
		}
	}
	return -1;
}

int main()
{
	while(scanf("%lld%lld%lld",&p,&a,&b)!=EOF)
	{
		ans=bsgs(a,b);
		if(ans==-1)printf("no solution\n");
		else printf("%lld\n",ans);
	}
	return 0;
}

• 例题：
  1.UVA10225 Discrete Logging（模板题）；
  2.【SDOI2013】随机数生成器；
  $$a^{i-1}\equiv \frac{x_i+\frac{b}{a+1}}{x_1+\frac{}{a-1}}$$
  其中$x_i=t$;
  3.【SDOI2011】计算器；

---

拓展大步小步算法

与普通的bsgs类似，但可以处理模数不为质数的情况。直接贴代码了：

#include<cstdio>
#include<cmath>
#include<map>
#define ll long long
using namespace std;

ll a,b,p,ans;
map<ll,ll>mp;
map<ll,bool>vis;

ll gcd(ll x,ll y)
{
	return (y==0)?x:gcd(y,x%y);
}

ll quick(ll x,ll y,ll mod)
{
	ll ret=1;
	while(y)
	{
		if(y&1)ret=ret*x%mod;
		x=x*x%mod;
		y>>=1;
	}
	return ret;
}

ll exbsgs(ll a,ll b,ll mod)
{
	if(b==1)return 0;
	ll k=0,tmp=1,d;
	while(1)
	{
		d=gcd(a,mod);
		if(d==1)break;
		if(b%d)return -1;
		b/=d;mod/=d;
		tmp=tmp*(a/d)%mod;
		k++;
		if(tmp==b)return k;
	}
	mp.clear();
	vis.clear();
	ll mul=b,q=ceil(sqrt(mod));
	mp[b]=0;
	for(int i=1;i<=q;++i)
	{
		mul=mul*a%mod;
		mp[mul]=i;
		vis[mul]=true;
	}
	ll qq=quick(a,q,mod);
	mul=tmp;
	for(int i=1;i<=q;++i)
	{
		mul=mul*qq%mod;
		if(vis[mul])return i*q-mp[mul]+k;
	}
	return -1;
}

int main()
{
	while(scanf("%lld%lld%lld",&p,&a,&b)!=EOF)
	{
		ans=exbsgs(a,b,p);
		if(ans==-1)printf("no solution\n");
		else printf("%lld\n",ans);
	}
	return 0;
}

---

快速乘和快速幂

• 应用：快速地做乘法或幂运算，多为辅助算法。
• 代码：

1.快速乘（乘法变加法）

#include<bits/stdc++.h>
using namespace std;
long long a,b,s,sum;
long long fast(long long a,long long b,long long mod)
{
	long long ans=0;
	while(b)
	{
		if(b%2)
		{
	    	ans+=a;	
	    	ans%=mod;
		}
		a+=a;
		a%=mod;
		b/=2;
	}
	return ans;
}
int main()
{
	scanf("%I64d%I64d%I64d",&a,&b,&s);
	sum=fast(a,b,s);
	printf("%I64d",sum);
	return 0;
}

2.快速幂：

#include<bits/stdc++.h>
using namespace std;
int  a,b,s,sum;
int fast(int a,int b,int mod)
{
	int ans=1;
	while(b)
	{
		if(b%2)
		{
	    	ans*=a;	
	    	ans%=mod;
		}
		a*=a;
		a%=mod;
		b=b>>1;
	}
	return ans;
}
int main()
{
	scanf("%d%d%d",&a,&b,&s);
	sum=fast(a,b,s);
	printf("%d",sum);
	return 0;
}

---

矩阵相关

• 矩阵
  这个就叫矩阵：
  $$\left[\begin{array}{r}{a}_{11},a_{12},...,a_{1n}\\ a_{21},a_{22},...,a_{2n}\\ ...,...,...,...\\ a_{n1},a_{n2},...,a_{nn}\end{array}\right]$$
• 矩阵运算

1.矩阵加减：
条件：两个矩阵长宽一样，例如
$$\left[\begin{array}{r}142\\ 200\end{array}\right]+\left[\begin{array}{r}005\\ 750\end{array}\right]=\left[\begin{array}{r}1+04+02+5\\ 2+70+50+0\end{array}\right]=\left[\begin{array}{r}147\\ 950\end{array}\right]$$
减法类似；
2.矩阵乘法：
条件：第一个矩阵的行数和第二个矩阵的列数相等。
$$\left[\begin{array}{r}102\\ -131\end{array}\right]\ast \left[\begin{array}{r}31\\ 21\\ 10\end{array}\right]=\left[\begin{array}{r}(1\ast 3+0\ast 2+2\ast 1)(1\ast 1+0\ast 1+2\ast 0)\\ (-1\ast 3+3\ast 2+1\ast 1)(-1\ast 1+3\ast 1+1\ast 0)\end{array}\right]=\left[\begin{array}{r}51\\ 42\end{array}\right]$$
3.矩阵的幂:
$${\left[\begin{array}{r}{a}_{11},a_{12},...,a_{1n}\\ a_{21},a_{22},...,a_{2n}\\ ...,...,...,...\\ a_{n1},a_{n2},...,a_{nn}\end{array}\right]}^2=\left[\begin{array}{r}{a}_{11},a_{12},...,a_{1n}\\ a_{21},a_{22},...,a_{2n}\\ ...,...,...,...\\ a_{n1},a_{n2},...,a_{nn}\end{array}\right]\ast \left[\begin{array}{r}{a}_{11},a_{12},...,a_{1n}\\ a_{21},a_{22},...,a_{2n}\\ ...,...,...,...\\ a_{n1},a_{n2},...,a_{nn}\end{array}\right]$$

• 矩阵快速幂
  是一种求矩阵的幂的方法，快速幂套矩阵乘法即可

#include<cstdio>
#include<iostream>
#define mod 1000000007
#define ll long long
using namespace std;
struct node
{
	ll m[105][105];
}a,e;//a是输入的矩阵，ans是答案矩阵，e是单位矩阵
ll n,k;
node mul(node x,node y)
{
	node nw;
	for(int i=1;i<=n;i++)
	for(int j=1;j<=n;j++)
	nw.m[i][j]=0;
	for(int i=1;i<=n;i++)
	for(int j=1;j<=n;j++)
	for(int k=1;k<=n;k++)
	{
		nw.m[i][j]=nw.m[i][j]%mod+x.m[i][k]*y.m[k][j]%mod;
	}
	return nw;
}//矩阵乘法
node quickk(node x,ll y)
{
	node pre=e;
	while(y)
	{
		if(y&1)
		pre=mul(pre,x);
		x=mul(x,x);
		y>>=1;
	}
	return pre;
}//快速幂主体
int main()
{
	scanf("%lld%lld",&n,&k);
	for(int i=1;i<=n;i++)
	for(int j=1;j<=n;j++)
	scanf("%lld",&a.m[i][j]);
	for(int i=1;i<=n;i++)
	e.m[i][i]=1;//构造单位矩阵
	node ans=quickk(a,k);
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
    	printf("%lld ",ans.m[i][j]%mod);
    	printf("\n");
	}
	return 0;
}

• 矩阵加速
  难点：构造矩阵
  若把运算放在矩阵中，会大大提高运算效率。
  举个例子：
  $a[1]=a[2]=a[3]=1$
  $a[x]=a[x-3]+a[x-1](x>3)$
  求$a$数列的第$n$项对1000000007$（10^9+7）$取余的值。
  $n$很大，如果递推求会超时。
  但是如果我们把这个模型放在矩阵中，求第二项就会变成求这个东西：
  $${\left[\begin{array}{r}100\\ 010\\ 101\end{array}\right]}^2$$
  求第$n$项：
  $${\left[\begin{array}{r}100\\ 010\\ 101\end{array}\right]}^n$$
  输出$(1,1)$。
  代码：

#include<cstdio>
#define ll long long
#define mod 1000000007
using namespace std;
struct node
{
	ll m[4][13];
}stp,e,ans;
ll n;
int t;
node mul(node x,node y)
{
	node nw;
	for(int i=1;i<=3;i++)
	for(int j=1;j<=3;j++)
	nw.m[i][j]=0;
	for(int i=1;i<=3;i++)
	for(int j=1;j<=3;j++)
	for(int k=1;k<=3;k++)
	{
		nw.m[i][j]=(nw.m[i][j]%mod+x.m[i][k]*y.m[k][j]%mod)%mod;
	}
	return nw;
}//矩阵乘法
node quickk(node x,ll y)
{
	node pre=e;
	while(y)
	{
		if(y&1)
		pre=mul(pre,x);
		x=mul(x,x);
		y>>=1;
	}
	return pre;
}//快速幂主体
int main()
{ 
    stp.m[1][14]=1;
    stp.m[1][15]=1;
    stp.m[2][16]=1;
    stp.m[3][17]=1;
	for(int i=1;i<=3;i++)
	e.m[i][i]=1;
	scanf("%d",&t);
	while(t--)
	{
		for(int i=1;i<=3;i++)
		for(int j=1;j<=3;j++)
		ans.m[i][j]=0;
    	scanf("%lld",&n);
    	ans=quickk(stp,n-1);
    	printf("%lld\n",ans.m[1][18]);
	}
	return 0;
}

例题：
1.【模板】矩阵加速（数列）
2.斐波那契数列

---

欧拉函数

• 概念：1~$x$-1中与$x$互质的数的个数。
• 特点：欧拉函数是积性函数。
• 代码：

1.线性求欧拉函数（求素数时顺便求）：

#include<cstdio>
#define maxn 1000005
using namespace std;

int n;
int phi[maxn],prime[maxn],tot,ans;
bool mark[maxn];

void getphi()
{
	phi[1]=1;
	for(int i=2;i<=n;++i)
	{
		if(!mark[i])
		{
			prime[++tot]=i;
			phi[i]=i-1;
		}//质数的phi
		for(int j=1;j<=tot;++j)
		{
			if(i*prime[j]>n)break;
			mark[i*prime[j]]=1;//标记合数
			if(i%prime[j]==0)
			{
				phi[i*prime[j]]=phi[i]*prime[j];
				break;
			}
			else phi[i*prime[j]]=phi[i]*(prime[j]-1);
		}
	}
}

int main()
{
	scanf("%d",&n);
	getphi();
	for(int i=1;i<=n;++i)
	printf("%d ",phi[i]);
	return 0;
}

2.非线性：

#include<bits/stdc++.h>
using namespace std;
int phi(int x)
{
	int ans=x;
	for(int i=2;i*i<=x;i++)
	{
		if(x%i==0)
		{
			ans=ans/i*(i-1);
			while(x%i==0)x/=i;
		}
	}
	if(x>1)ans=ans/x*(x-1);//x为质数
	return ans;
}
int main()
{
	int n;
	while(scanf("%d",&n)==1)
	printf("%d\n",phi(n));
	return 0;
}

---

欧拉定理及费马小定理

• 欧拉定理：若$a$与$m$互质，则$a^{\phi (m)}\equiv 1(modm)$;
• 拓展欧拉定理：若$r>\phi (m)$，则$a^r\equiv a^{rmod\phi (m)+\phi (m)}(modm)$;
• 费马小定理：$a^{p-1}\equiv 1(modq)$，其中$q$为质数。（费马小定理是欧拉定理的特殊形式）
• 代码
  欧拉定理：

#include<cstdio>
#define ll long long
using namespace std;

ll a,b,m;
ll tmp,phi,flag;
char ch;

ll quick(ll x,ll y)
{
	ll ans=1;
	while(y)
	{
		if(y&1)ans=ans*x%m;
		x=x*x%m;
		y>>=1;
	}
	return ans;
}

int main()
{
	scanf("%lld%lld",&a,&m);
	tmp=phi=m;
	for(ll i=2;i*i<=m;++i)
	{
		if(tmp%i)continue;
		phi-=phi/i;
		while(tmp%i==0)
		tmp/=i;
	}
	if(tmp>1)
	phi-=phi/tmp;
	while((ch=getchar())<'0'||ch>'9');
	b=ch-'0';
	if(b>=phi)
	{
		flag=1;
		b%=phi;
	}
	while((ch=getchar())>='0'&&ch<='9')
	{
		b=b*10+ch-'0';
		if(b>=phi)
		{
			flag=1;
			b%=phi;
		}
	}
	if(flag)b+=phi;
	printf("%lld",quick(a,b));
	return 0;
}

模板题：【模板】欧拉定理

---

中国剩余定理

• 作用：用于求解同余方程组：
  $$\{\begin{array}{r}x\equiv b_1(mod{m}_1)\\ x\equiv b_2(mod{m}_2)\\ ................\\ x\equiv b_n(mod{m}_n)\end{array}$$
  其中$m_1,m_2,...,m_n$互质。
• 代码

#include<cstdio>
#define maxn 100005
using namespace std;

int n;
int a[maxn],m[maxn];

void exgcd(int a,int b,int &x,int &y)
{
	if(b==0)
	{
		x=1;
		y=0;
		return;
	}
	int c=a-(a/b)*b,xx,yy;
	exgcd(b,c,xx,yy);
	x=yy;y=xx-(a/b)*yy;
}

int china(int w[],int b[],int k)
{
	int x,y,a=0,m,n=1;
	for(int i=1;i<=k;++i)
	n*=w[i];
	for(int i=1;i<=k;++i)
	{
		m=n/w[i];
		exgcd(w[i],m,x,y);
		a=(a+y*m*b[i])%n;
	}
	if(a>0)return a;
	else return a+n;
}

int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;++i)
	scanf("%d%d",&m[i],&a[i]);
	printf("%d",china(m,a,n));
	return 0;
}

---

拓展中国剩余定理

• 作用
  用于解同余方程组，但模数两两之间可以不互质。
• 证明
  我写的博客：拓展中国剩余定理
• 代码

#include<cstdio>
#define ll long long
using namespace std;
ll n,a,b,tot,M;
ll ans,x,y,r;
ll flag=1,mul;
ll exgcd(ll a,ll b,ll &x,ll &y)
{
	if(b==0)
	{
		x=1;y=0;
		return a;
	}
	ll d=exgcd(b,a%b,x,y);
	ll t=y;
	y=x-(a/b)*y;
	x=t;
	return d;
}//求解同余方程
ll quick_mul(ll nw,ll aim,ll mod)
{
	ll res=0;
	while(aim>0)
	{
		if(aim&1)res=(res+nw)%mod;
		nw=(nw+nw)%mod;
		aim>>=1;
	}
	return res;
}//龟速乘QAQ（模板）
int main()
{
	scanf("%I64d",&n);
	while(n--)
	{
		x=0;
		y=0;
		tot++;//tot计数
		scanf("%I64d%I64d",&a,&b);
		if(tot==1)
		{
			ans=b;
			M=a;
			continue;
		}//初始化
		r=exgcd(M,a,x,y);//求最大公因数
		mul=((b-ans)%a+a)%a;
		x=quick_mul(x,mul/r,a);
		if((b-ans)%r!=0)
		{
			flag=0;
			continue;
		}//判断无解
		ans=ans+(x*M);
		M=(M*a)/r;//更新最小公倍数
		ans=(ans%M+M)%M;//保证ans为正值
	}
	if(flag)
	printf("%I64d",ans);
	else
	printf("No");
	return 0;
}

---

卢卡斯定理

• 作用
  求$C_n^m\%p$($p$为质数)的值。
• 定理：
  1.$$C_n^m=C_{n/p}^{m/p}\ast C_{n\%p}^{m\%p}$$
  2.将$n,m$写成$p$进制：
  $$C_n^m=C_{n_{k-1}}^{m_{k-1}}\ast C_{n_{k-2}}^{m_{k-2}}\ast ...\ast C_{n_0}^{m_0}(modp)$$
• 代码：

#include<cstdio>
#define ll long long
using namespace std;

ll t,n,m,p;

ll quick(ll x,ll y)
{
	ll ret=1;
	while(y)
	{
		if(y&1)ret=(ret*x)%p;
		x=(x*x)%p;
		y>>=1;
	}
	return ret;
}

ll C(ll n,ll m)
{
	if(n<m)return 0;
	if(m>n-m)m=n-m;
	ll s1=1,s2=1;
	for(ll i=0;i<m;++i)
	{
		s1=s1*(n-i)%p;
		s2=s2*(i+1)%p;
	}
	return s1*quick(s2,p-2)%p;
}

ll lucas(ll n,ll m)
{
	if(m==0)return 1;
	return C(n%p,m%p)*lucas(n/p,m/p)%p;
}

int main()
{
	scanf("%lld",&t);
	while(t--)
	{
		scanf("%lld%lld%lld",&n,&m,&p);
		printf("%lld\n",lucas(m+n,m));
	}
	return 0;
}

模板题：【模板】卢卡斯定理

---

拓展卢卡斯定理

• 作用：同卢卡斯定理，只是$p$不一定为质数。但其实与卢卡斯定理没什么关系。
• 代码（我证不来$QAQ$）:

#include<cstdio>
#include<cmath>
#define ll long long
#define maxn 1000005
using namespace std;

ll n,m,p;
ll A[maxn],B[maxn];

ll gcd(ll a,ll b)
{
	return (b==0)?a:gcd(b,a%b);
}

void exgcd(ll a,ll b,ll &x,ll &y)
{
	if(b==0)
	{
		x=1;
		y=0;
		return;
	}
	ll c=a-(a/b)*b,xx,yy;
	exgcd(b,c,xx,yy);
	x=yy;y=xx-(a/b)*yy;
}

ll inv(ll a,ll b)
{
	ll x,y;
	exgcd(a,b,x,y);
	return (x+b)%b;
} 

ll quick(ll x,ll y,ll mod)
{
	ll ret=1;
	while(y)
	{
		if(y&1)ret=(ret*x)%mod;
		x=(x*x)%mod;
		y>>=1;
	}
	return ret;
}

ll f(ll n,ll x,ll y)
{
	if(n==0)return 1;
	ll ret=1;
	for(ll i=2;i<=y;++i)
	{
		if(i%x)
		ret=(ret*i)%y;
	}
	ret=quick(ret,n/y,y);
	for(ll i=2;i<=n%y;++i)
	{
		if(i%x)
		ret=(ret*i)%y;
	}
	return ret*f(n/x,x,y)%y;
}

ll china(ll x,ll y)
{
	return x*inv(p/y,y)%p*(p/y)%p;
}

ll C(ll n,ll m,ll x,ll y)
{
	ll f1=f(n,x,y),f2=f(m,x,y),f3=f(n-m,x,y);
	ll ret=0;
	for(ll i=n;i;i/=x)ret+=i/x;
	for(ll i=m;i;i/=x)ret-=i/x;
	for(ll i=n-m;i;i/=x)ret-=i/x;
	return f1*inv(f2,y)%y*inv(f3,y)%y*quick(x,ret,y)%y;
}

ll exlucas(ll n,ll m,ll mod)
{
	ll ret=0,q=sqrt(mod)+5,x,tmp=mod;
	for(ll i=2;i<=q;++i)
	{
		if(tmp%i)continue;
		x=1;
		while(tmp%i==0)
		{
			x*=i;
			tmp/=i;
		}
		ret=(ret+china(C(n,m,i,x),x))%mod;
	}
	if(tmp>1)ret=(ret+china(C(n,m,tmp,tmp),tmp))%mod;
	return ret;
}

int main()
{
	scanf("%lld%lld%lld",&n,&m,&p);
	printf("%lld",exlucas(n,m,p));
	return 0;
}

---

狄利克雷卷积

• 概念
  定义一种运算:
  $$h=f⨂g$$
  为：
  $$h(n)=\sum _{d|n}f(d)g(\frac{d}{n})$$
  这就是狄利克雷卷积的一种表现形式（详情见博客：Greninja_Wu 狄利克雷卷积）。
  若$f$,$g$为积性函数，则$h$也为积性函数。

• 应用
  莫比乌斯反演。
  本章完。

---

莫比乌斯函数

• 概念
  定义：
  $$\mu (n)=\{\begin{array}{r}0,p^2|n\\ (-1{)}^r,n=p_1\ast p_2\ast ...\ast p_r\end{array}$$
• 作用
  常用于与约数相关的容斥。
• 求法
  可以线性筛预处理（同欧拉函数）。
• 代码

#include<cstdio>
#define maxn 1000005
using namespace std;

int n;
int mu[maxn],prime[maxn],tot,ans;
bool mark[maxn];

void getmu()
{
	mu[1]=1;
	for(int i=2;i<=n;++i)
	{
		if(!mark[i])
		{
			prime[++tot]=i;
			mu[i]=-1;
		}//质数的phi
		for(int j=1;j<=tot;++j)
		{
			if(i*prime[j]>n)break;
			mark[i*prime[j]]=1;//标记合数
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				break;
			}
			else mu[i*prime[j]]=-mu[i];
		}
	}
}

int main()
{
	scanf("%d",&n);
	getmu();
	for(int i=1;i<=n;++i)
	printf("%d ",mu[i]);
	return 0;
}

---

莫比乌斯反演

• 前置定理：
  1.$$\sum _{d|n}\mu (d)=[n=1]$$
  也就是说，当$n=1$时，上式的值为1，否则为0.
  2.$$\sum _{d|n}\phi (d)=n$$
• 结论
  给定数论函数$f$和$g$：
  若$f$满足$f(n)={\sum }_{d|n}g(d)$，则$g$满足$ g(n)=\sum_{d|n}\mu(d)*f(\frac{n}{d})$
• 证明
  ${\sum }_{d|n}\mu (d)\ast f(\frac{n}{d})={\sum }_{d|n}\mu (d){\sum }_{d^{{}^{\prime }}|\frac{n}{d}}g(d^{{}^{\prime }})={\sum }_{d^{{}^{\prime }}|n}g(d^{{}^{\prime }}){\sum }_{d|\frac{n}{d^{{}^{\prime }}}}\mu (d)=g(n)$
• 应用
  将很复杂的式子变形成为很简单的式子。
• 例题
  1.【SDOI2015】约数个数和；
  2.【POI2007】ZAP-Queries；
  双倍经验：
  3.GCD；
  4.YY的GCD；
• 总结
  这种题通常是先把复杂的式子化为简单的式子，再用数论分块来做。
  这里以ZAP为例：
  目标公式：$$\sum _{i=1}^m\sum _{j=1}^n[gcd(i,j)=k]$$
  转化：
  ${\sum }_{i=1}^m{\sum }_{j=1}^n[gcd(i,j)=k]$
  $={\sum }_{i=1}^m{\sum }_{j=1}^n{\sum }_{d|gcd(i,k)/k}\mu (d)$
  $={\sum }_d\mu (d){\sum }_{kd|i}^m{\sum }_{kd|j}^{n}1$
  $={\sum }_d\mu (d)\lfloor \frac{m}{kd}\rfloor \lfloor \frac{n}{kd}\rfloor$;
  预处理$\lfloor \frac{m}{kd}\rfloor$和$\lfloor \frac{n}{kd}\rfloor$
• 代码：

#include<cstdio>
#include<algorithm>
#define ll long long
#define maxn 50005
using namespace std;

ll t,a,b,d,tot;
ll prime[maxn],mu[maxn],sum[maxn];
bool mark[maxn];

void getmu()
{
	mu[1]=1;
	for(ll i=2;i<=50000;++i)
	{
		if(!mark[i])
		{
			prime[++tot]=i;
			mu[i]=-1;
		}
		for(ll j=1;j<=tot;++j)
		{
			if(i*prime[j]>50000)break;
			mark[i*prime[j]]=1;
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				break;
			}
			else mu[i*prime[j]]=-mu[i];
		}
	}
	for(ll i=1;i<=50000;++i)
	sum[i]=sum[i-1]+mu[i];
}

int main()
{
	getmu();
	scanf("%lld",&t);
	while(t--)
	{
		ll n,ans=0;
		scanf("%lld%lld%lld",&a,&b,&d);
		n=min(a,b);
		for(ll l=1,r;l<=n;l=r+1)
		{
			r=min(a/(a/l),b/(b/l));
			ans+=(a/(l*d))*(b/(l*d))*(sum[r]-sum[l-1]);
		}
		printf("%lld\n",ans);
	}
	return 0;
}

---

杜教筛

• 应用
  求前缀和。
  其实求前缀和可以用线性的方法，但是当$n$很大（如$10^{10}$）时，就必须使用杜教筛了。
• 实质
  记忆化搜索。
• 实现方法
  定义数论函数$f$,$g$,$h$，其中$h=f⨂g$。
  目标是求$f$的前缀和$S(n)={\sum }_{i=1}^{n}f(i)$，而$h$的前缀和易求，可用此方法。
• 关键：找$g$和$h$
• 过程
  ${\sum }_{i=1}^{n}h(i)$
  $={\sum }_{i=1}^n{\sum }_{j|i}f(j)g(\frac{i}{j})$(狄利克雷卷积)
  $={\sum }_{i=1}^{n}g(i){\sum }_{j=1}^{\lfloor n/i\rfloor }f(j)$（玄学变换）
  $={\sum }_{i=1}^{n}g(i)S(\lfloor \frac{n}{i}\rfloor )$
  $=g(1)S(n)+{\sum }_{i=2}^{n}g(i)S(\lfloor \frac{n}{i}\rfloor )$

即：
$$S(n)=\frac{{\sum }_{i=1}^{n}h(i)-{\sum }_{i=2}^{n}g(i)S(\lfloor \frac{n}{i}\rfloor )}{g(1)}$$
用数论分块即可求出正解。
总复杂度:$O(n^{\frac{3}{4}})$

• 模板
  先预处理出小范围的答案，大范围的用map存。
  求${\sum }_{i=1}^n\mu (i)$

#include<cstdio>
#include<algorithm>
#include<map>
#define ll long long
#define maxn 5000005
#define mod 1000000007
using namespace std;

map<ll,ll>mp;
ll n,m,tot;
ll ans[maxn],mu[maxn],prime[maxn];
bool mark[maxn];

void getmu()
{
	mu[1]=1;
	for(ll i=2;i<=m;++i)
	{
		if(!mark[i])
		{
			prime[++tot]=i;
			mu[i]=-1;
		}//质数的phi
		for(ll j=1;j<=tot;++j)
		{
			if(i*prime[j]>m)break;
			mark[i*prime[j]]=1;//标记合数
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				break;
			}
			else mu[i*prime[j]]=-mu[i];
		}
	}
}

inline ll f(ll n)
{
	if(n<=5000000)return ans[n];
	if(mp.find(n)!=mp.end())return mp[n];
	ll ret=1;
	for(ll l=2,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		ret-=(r-l+1)*f(n/l);
		ret=(ret%mod+mod)%mod;
	}
	return mp[n]=ret%mod;
}

int main()
{
	scanf("%lld",&n);
	m=min(n,(ll)5000000);
	getmu();
	for(ll i=1;i<=m;++i)
	{
		ans[i]=(ans[i-1]+i*mu[i])%mod;
		ans[i]=(ans[i]+mod)%mod;
	}
	printf("%lld",f(n));
	return 0;
}

• 例题
  【模板】杜教筛（卡常好题）。

---

快速傅里叶变换（FFT）

• 复数
  1.形如：$a+bi$的数，其中$i^2=-1$。
  2.四则运算：
  a.$(a+bi)+(c+di)=(a+c)+(b+d)i$;
  b.$(a+bi)-(c+di)=(a-c)+(b-d)i$;
  c.$(a+bi)\ast (c+di)=(ac-bd)+(bc+ad)i$;
  d.除法被吞了$qwq$;

• 单位复数根：
  1.欧拉公式:$e^{ix}=\cos x+i\sin x$；
  2.设${\omega }_n^k$表示$n$的$k$次复根，即把一个平面均分成$n$份，取逆时针前$k$份。当$n=8$时，${\omega }_n^k$如下：
  
  3.根据欧拉公式：${\omega }_n^k=\cos (\frac{2\pi k}{n})+i\sin (\frac{2\pi k}{n})$
  4.性质：
  a.${\omega }_n^n={\omega }_n^0=1$;
  b.消去引理：${\omega }_{dn}^{dk}={\omega }_n^k$;
  c.折半引理：$({\omega }_n^k{)}^2={\omega }_{\frac{n}{2}}^k$;
  d.${\omega }_{n_1}^{k_1}\ast {\omega }_{n_2}^{k_2}={\omega }_{n_1+n_2}^{k_1+k_2}$;

• 快速傅里叶变换
  终于进入正题了，写的我好辛苦啊！！！
  1.应用：用于求多项式的乘积，甚至可以求$A\ast BProblem$(雾)。
  2.算法：Cooley-Tukey（分治）
  一般求多项式的乘积用的是系数表示法$A={\sum }_{i=0}^n{a}_i\ast x^i$，这样求多项式的乘积的时间复杂度就是$O(n^2)$，有点慢。
  而FFT选用的是点值表示法$(x_1,y_1),...,(x_n,y_n)$，多项式的乘积就是$(x_1,y_{a_1}\ast y_{b_1}),...,(x_n,y_{a_n}\ast y_{b_n})$，时间复杂度是$O(n)$，好快啊！！！
  接下来就是选择点值表达式中的点了。上面我们提到了单位复数根，它有很多优美的性质，我们就选它作为点横坐标啦(雾)。
  那么:$$A({\omega }_n^k)=\sum _{i=0}^{n-1}{a}_i\ast {\omega }_n^{k_i}$$
  把$A(x)$的系数按奇偶分开，即：
  $$A_0(x)=a_0+a_{2}x+a_4{x}^2+...$$
  $$A_1(x)=a_1+a_{3}x+a_5{x}^2+...$$
  $$A(x)=A_0(x^2)+x{A}_1(x^2)$$
  有兴趣可以自己验证一下。
  那么：$$A({\omega }_n^k)=\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i}{\omega }_n^{2ki}+{\omega }_n^k\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{\omega }_n^{2ki}$$
  当$k<\frac{n}{2}$时：
  $$A({\omega }_n^k)=\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i}{\omega }_{\frac{n}{2}}^{ki}+{\omega }_n^k\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{\omega }_{\frac{n}{2}}^{ki}$$
  $$A({\omega }_n^{k+\frac{n}{2}})=\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i}{\omega }_{\frac{n}{2}}^{ki}+{\omega }_n^k\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{\omega }_{\frac{n}{2}}^{ki}$$
  每次代入规模减半，复杂度为$O(n\log n)$

• 傅里叶逆变换
  将点值表达式转化为系数表达式。

• 实现
  1.递归实现（常数太大，不常用）；
  2.非递归实现：先位翻转(蝴蝶变换)，再自底向上实现，常数较小。

• 代码
  其实说了那么多，我也不知道自己在写什么，还是看代码吧。

#include<cstdio>
#include<cmath>
#include<algorithm>
#define maxn 10000005
using namespace std;

const double pi=acos(-1.0);
struct complex
{
	double x,y;
	complex(double xx=0,double yy=0){x=xx,y=yy;}
}a[maxn],b[maxn];
int n,m;
int lim=1,l,r[maxn];

complex operator+(complex a,complex b)
{
	return complex(a.x+b.x,a.y+b.y);
}

complex operator-(complex a,complex b)
{
	return complex(a.x-b.x,a.y-b.y);
}

complex operator*(complex a,complex b)
{
	return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}//复数运算

void fft(complex *c,int flag)
{
	for(int i=0;i<lim;++i)
	{
		if(i<r[i])
		swap(c[i],c[r[i]]);
	}//求出迭代序列
	for(int mid=1;mid<lim;mid<<=1)
	{
		complex wn(cos(pi/mid),flag*sin(pi/mid));
		for(int r=mid<<1,j=0;j<lim;j+=r)//r是右端点,j表示当前位置
		{
			complex w(1,0);
			for(int k=0;k<mid;++k,w=w*wn)
			{
				complex x=c[j+k],y=w*c[j+mid+k];
				c[j+k]=x+y;
				c[j+mid+k]=x-y;
			}
		}
	}
}

int main()
{
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;++i)
	scanf("%lf",&a[i].x);
	for(int i=0;i<=m;++i)
	scanf("%lf",&b[i].x);
	while(lim<=n+m)
	{
		lim<<=1;
		l++;
	}
	for(int i=0;i<lim;++i)
	r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));//蝴蝶变换
	fft(a,1);
	fft(b,1);
	for(int i=0;i<=lim;++i)
	a[i]=a[i]*b[i];
	fft(a,-1);//傅里叶逆变换
	for(int i=0;i<=n+m;++i)
	printf("%d ",(int)(a[i].x/lim+0.5));
	return 0;
}

• 总结
  形如求$c(n)={\sum }_{i=1}^{n-1}a[i]+b[n-i]$的式子适合用$fft$或$ntt$解决。

• 例题
  1.【模板】多项式乘法（FFT）；
  2.【模板】分治 FFT；
  3.【模板】A*B Problem升级版（FFT快速傅里叶）；

---

快速数论变换（NTT）

• 原根
  若$g$是$p$的一个原根，则$g^0,g^1,...,g^{p-1}(modp)$是$p-1$个非0数。

• NTT
  与FFT类似，只是有模，并且${\omega }_n^i$变成了原根$g$而已。

• 代码
  改天再贴。

---

*斯特林数相关算法

由于太难了，有时间再写。