算法笔记——数学相关

算法 数学

---

高精度

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

• 应用范围：

  1.$long long$存不下$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*\frac{1}{x}=1$，在$c++$中如果贸然直接除以$x$，可能会出问题。然而如果对于一个数$x$，若存在另一个数$x^{-1}$，使$x*x^{-1}$模上一个数$mod$等于1，那么这个数$x^{-1}$就是$x$的逆元。

• 求法

1.递推法：$inv[i]=-(p/i)*inv[p mod i]$(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^n_m=(m-n)!$

2.错位排列（错排）
从$m$个数选出$n$个数，数$n$不能在第$n$个位置上的方案数。
$f(n)=(n-1)*(f(n-1)+f(n-2))$
3.圆排列
从$n$个数选出$r$个数围成一个圈的方案数。

$A^r_n/r$,若$r=n$，则有(n-1)!种。

• 组合

$C^n_m=\frac{m!}{n!*(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\leq m}_{i=0}(^i_m)x^iy^{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*prime[j]]$标记为合数，若$i$%$prime[j]$=0，说明$i$是$prime[j]$的倍数，跳出循环，保证每一个合数只会被筛一次。
举个例子，当$i$为4时，质数数组有2和3，枚举质数2，标记4*2=8为合数。因为4是2的倍数，跳出循环。假如不跳出循环，那么下一个被枚举到的合数就是4*3=12，而12应该在枚举6*2时筛到，多筛了一遍，增加了时间复杂度。

4.Miller-Rabin素性测试：
玄学算法。
费马小定理：若$p$是质数，则$\forall x^{p-1}\equiv1(mod p)$，但不是所有数只要满足$\forall x^{p-1}\equiv1(mod p)$就是质数，所以光有这一个定理还不够;
二次探测定理：若$p$是质数，且$\forall x^2\equiv1(mod p)$，则$(x+1)(x-1)\equiv0(mod p)$，也就是说$x\equiv\pm1(mod p)$。
Miller-Rabin的流程：
设$p$是大于2的奇数，另$p-1=d*2^r$，d为奇数；
随机选择一个$x$，
首先判断费马小定理是否成立：$x^{p-1}\equiv =1(mod p)$；
若符合，则由二项式定理，要么$x^d\equiv1(mod p)$，要么存在一个比$r$小的数$i$满足$x^{d*2^i}\equiv-1(mod p)$。
单次复杂度为$O(log n)$，多次随机$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^{r_i}_{k=0}p_i^{}$;
4.约数函数：$\sigma_t(n)=\sum_{d|n}dt=\prod \sum^{r_i}_{k=0}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(mod b)$，以及不定方程$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*a+0*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;
}

其实，求解同余方程，还有一种方法：

a%=p;b%=p;
if(a==0&&b!=0)
printf("Orz, I cannot find x!\n");
else
printf("%lld\n",quick(a,p-2,p)*b%p);
continue;

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

• 例题：

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

---

大步小步算法（BSGS）

• 应用：求解高次同余方程
  $$a^x \equiv b(mod\ c)$$
  其中$c$为质数。
• 由于我也证不来为什么算法是对的，所以我直接粘贴代码了$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{aligned} a_{11},a_{12},...,a_{1n}\\ a_{21},a_{22},...,a_{2n}\\ ...,...,...,...\\ a_{n1},a_{n2},...,a_{nn} \end{aligned}\right]$$
• 矩阵运算

1.矩阵加减：
条件：两个矩阵长宽一样，例如
$\left[\begin{aligned}1\ 4\ 2\\2\ 0\ 0\end{aligned}\right]+\left[\begin{aligned}0\ 0\ 5\\7\ 5\ 0\end{aligned}\right]=\left[\begin{aligned}1+0\ 4+0\ 2+5\\2+7\ 0+5\ 0+0\end{aligned}\right]=\left[\begin{aligned}1\ 4\ 7\\9\ 5\ 0\end{aligned}\right]$
减法类似；
2.矩阵乘法：
条件：第一个矩阵的行数和第二个矩阵的列数相等。
$\left[\begin{aligned}1\ 0\ 2\\-1\ 3\ 1\end{aligned}\right]*\left[\begin{aligned}3\ 1\\ 2\ 1\\1\ 0\end{aligned}\right]=\left[\begin{aligned}(1*3+0*2+2*1)\ (1*1+0*1+2*0)\\ (-1*3+3*2+1*1)\ (-1*1+3*1+1*0)\end{aligned}\right]=\left[\begin{aligned}5\ 1\\ 4\ 2\end{aligned}\right]$
3.矩阵的幂:
$\left[\begin{aligned} a_{11},a_{12},...,a_{1n}\\ a_{21},a_{22},...,a_{2n}\\ ...,...,...,...\\ a_{n1},a_{n2},...,a_{nn} \end{aligned}\right]^2=\left[\begin{aligned} a_{11},a_{12},...,a_{1n}\\ a_{21},a_{22},...,a_{2n}\\ ...,...,...,...\\ a_{n1},a_{n2},...,a_{nn} \end{aligned}\right]*\left[\begin{aligned} a_{11},a_{12},...,a_{1n}\\ a_{21},a_{22},...,a_{2n}\\ ...,...,...,...\\ a_{n1},a_{n2},...,a_{nn} \end{aligned}\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{aligned}1\ 0\ 0\\0\ 1\ 0\\1\ 0\ 1 \end{aligned}\right]^2$$
  求第$n$项：
  $$\left[\begin{aligned}1\ 0\ 0\\0\ 1\ 0\\1\ 0\ 1 \end{aligned}\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^{\varphi(m)}\equiv 1(mod\ m)$;
• 拓展欧拉定理：若$r>\varphi(m)$，则$a^r\equiv a^{r\ mod\ \varphi(m)+\varphi(m)}(mod\ m)$;
• 费马小定理：$a^{p-1}\equiv 1(mod\ q)$，其中$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;
}

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

---

中国剩余定理

• 作用：用于求解同余方程组：
  $$\left\{\begin{aligned}x\equiv b_1(mod\ m_1)\\x\equiv b_2(mod\ m_2)\\................\\x\equiv b_n(mod\ m_n)\end{aligned}\right.$$
  其中$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^m_n\%p$($p$为质数)的值。
• 定理：
  1. $$C^m_n=C^{m/p}_{n/p}*C^{m\%p}_{n\%p}$$
  2.将$n,m$写成$p$进制：
  $$C^m_n=C^{m_{k-1}}_{n_{k-1}}*C^{m_{k-2}}_{n_{k-2}}*...*C^{m_0}_{n_0}(mod\ p)$$
• 代码：

#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\bigotimes g$$
  为：
  $$h(n)=\sum_{d|n}f(d)g(\frac{d}{n})$$
  这就是狄利克雷卷积的一种表现形式（详情见博客：Greninja_Wu 狄利克雷卷积）。
  若$f$,$g$为积性函数，则$h$也为积性函数。

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

---

莫比乌斯函数

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

#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}\varphi(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)*f(\frac{n}{d})=\sum_{d|n}\mu(d)\sum_{d^{'}|\frac{n}{d}}g(d^{'})=\sum_{d^{'}|n}g(d^{'})\sum_{d|\frac{n}{d^{'}}}\mu(d)=g(n)$
• 应用
  将很复杂的式子变形成为很简单的式子。
• 例题
  1.【SDOI2015】约数个数和；
  2.【POI2007】ZAP-Queries；
  双倍经验：
  3.GCD；
  4.YY的GCD；
• 总结
  这种题通常是先把复杂的式子化为简单的式子，再用数论分块来做。
  这里以ZAP为例：
  目标公式： $$\sum^m_{i=1}\sum^n_{j=1}[gcd(i,j)=k]$$
  转化：
  $\sum^m_{i=1}\sum^n_{j=1}[gcd(i,j)=k]$
  $=\sum^m_{i=1}\sum^n_{j=1}\sum_{d|gcd(i,k)/k}\mu(d)$
  $=\sum_d\mu(d)\sum^m_{kd|i}\sum^n_{kd|j}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\bigotimes g$。
  目标是求$f$的前缀和$S(n)=\sum^n_{i=1}f(i)$，而$h$的前缀和易求，可用此方法。
• 关键：找$g$和$h$
• 过程
  $\sum^n_{i=1}h(i)$
  $=\sum^n_{i=1}\sum_{j|i}f(j)g(\frac{i}{j})$(狄利克雷卷积)
  $=\sum^n_{i=1}g(i)\sum^{\lfloor n/i \rfloor}_{j=1}f(j)$（玄学变换）
  $=\sum^n_{i=1}g(i)S(\lfloor \frac{n}{i} \rfloor)$
  $=g(1)S(n)+\sum^n_{i=2}g(i)S(\lfloor \frac{n}{i} \rfloor)$

即：

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

• 模板
  先预处理出小范围的答案，大范围的用map存。
  求$\sum^n_{i=1} \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)*(c+di)=(ac-bd)+(bc+ad)i$;
  d.除法被吞了$qwq$;

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

• 快速傅里叶变换
  终于进入正题了，写的我好辛苦啊！！！
  1.应用：用于求多项式的乘积，甚至可以求$A*B\ Problem$(雾)。
  2.算法：Cooley-Tukey（分治）
  一般求多项式的乘积用的是系数表示法$A=\sum^n_{i=0}a_i*x^i$，这样求多项式的乘积的时间复杂度就是$O(n^2)$，有点慢。
  而FFT选用的是点值表示法$(x_1,y_1),...,(x_n,y_n)$，多项式的乘积就是$(x_1,y_{a_1}*y_{b_1}),...,(x_n,y_{a_n}*y_{b_n})$，时间复杂度是$O(n)$，好快啊！！！
  接下来就是选择点值表达式中的点了。上面我们提到了单位复数根，它有很多优美的性质，我们就选它作为点横坐标啦(雾)。
  那么:

  $$A(\omega^k_n)=\sum^{n-1}_{i=0}a_i*\omega^{k_i}_n$$
  把$A(x)$的系数按奇偶分开，即：
  $$A_0(x)=a_0+a_2x+a_4x^2+...$$
  $$A_1(x)=a_1+a_3x+a_5x^2+...$$
  $$A(x)=A_0(x^2)+xA_1(x^2)$$
  有兴趣可以自己验证一下。
  那么： $$A(\omega^k_n)=\sum^{\frac{n}{2}-1}_{i=0}a_{2i}\omega^{2ki}_{n}+\omega^k_n\sum^{\frac{n}{2}-1}_{i=0}a_{2i+1}\omega^{2ki}_{n}$$
  当$k<\frac{n}{2}$时：
  $$A(\omega^k_n)=\sum^{\frac{n}{2}-1}_{i=0}a_{2i}\omega^{ki}_{\frac{n}{2}}+\omega^k_n\sum^{\frac{n}{2}-1}_{i=0}a_{2i+1}\omega^{ki}_{\frac{n}{2}}$$
  $$A(\omega^{k+\frac{n}{2}}_n)=\sum^{\frac{n}{2}-1}_{i=0}a_{2i}\omega^{ki}_{\frac{n}{2}}+\omega^k_n\sum^{\frac{n}{2}-1}_{i=0}a_{2i+1}\omega^{ki}_{\frac{n}{2}}$$
  每次代入规模减半，复杂度为$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^{n-1}_{i=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}(mod\ p)$是$p-1$个非0数。

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

• 代码
  改天再贴。

---

*斯特林数相关算法

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