---

个人算法笔记：戳我qwq

高精度

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

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

• 求法

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

• 应用：

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;
}//快速幂求逆元

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

---

排列组合

• 排列

排列的分类：
1.选排列
最普通的排列，从mmm个数选出nnn个数的方案数。
Amn=(m−n)!A^n_m=(m-n)!Amn​=(m−n)!
2.错位排列（错排）
从mmm个数选出nnn个数，数nnn不能在第nnn个位置上的方案数。
f(n)=(n−1)∗(f(n−1)+f(n−2))f(n)=(n-1)*(f(n-1)+f(n-2))f(n)=(n−1)∗(f(n−1)+f(n−2))
3.圆排列
从nnn个数选出rrr个数围成一个圈的方案数。
Anr/rA^r_n/rAnr​/r,若r=nr=nr=n，则有(n-1)!种。

• 组合
  Cmn=m!n!∗(m−n)!C^n_m=\frac{m!}{n!*(m-n)!}Cmn​=n!∗(m−n)!m!​
  模板：逆元求组合数

#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=∑i=0i≤m(mi)xiym−i(x+y)^m=\sum^{i\leq m}_{i=0}(^i_m)x^iy^{m-i}(x+y)m=∑i=0i≤m​(mi​)xiym−i

---

质数的判定和应用

• 质数的判定方法

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

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

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

• 代码

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.将nnn分解质因数：n=∏pirin=\prod p_i^{r_i}n=∏piri​​;
2.约数个数：σ0(n)=∏(1+ri)\sigma_0(n)=\prod (1+r_i)σ0​(n)=∏(1+ri​);
3.约数和：σ1(n)=∏∑k=0ripi\sigma_1(n)=\prod \sum^{r_i}_{k=0}p_i^{}σ1​(n)=∏∑k=0ri​​pi​;
4.约数函数：σt(n)=∑d∣ndt=∏∑k=0ripikt\sigma_t(n)=\sum_{d|n}dt=\prod \sum^{r_i}_{k=0}p_i^{kt}σt​(n)=∑d∣n​dt=∏∑k=0ri​​pikt​;
以上这些东西，我也看不懂。
例题：
1.[AHOI2007]密码箱（数据过水）;
2.[HAOI2007]反素数;
3.[JLOI2014]聪明的燕姿;

• 最大公因数(gcd)

最大公因数的性质：
1.gcd(a,b)=gcd(a−b,b)gcd(a,b)=gcd(a-b,b)gcd(a,b)=gcd(a−b,b);
2.gcd(a,b)=gcd(a%b,b)gcd(a,b)=gcd(a\%b,b)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≡m(modb)ax\equiv m(mod b)ax≡m(modb)，以及不定方程ax+by=gcd(a,b)ax+by=gcd(a,b)ax+by=gcd(a,b)。

• 方法

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

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

• 步骤

斐蜀定理：ax+by=gcd(a,b)ax+by=gcd(a,b)ax+by=gcd(a,b)一定有解。
gcd(a,b)=gcd(b,a%b)gcd(a,b)=gcd(b,a\%b)gcd(a,b)=gcd(b,a%b)。
在欧几里得算法的最后一步时，有b=0b=0b=0，gcd(a,b)=agcd(a,b)=agcd(a,b)=a，因为1∗a+0∗0=a1*a+0*0=a1∗a+0∗0=a，所以此时ax+by=gcd(a,b)ax+by=gcd(a,b)ax+by=gcd(a,b)有一组解x=1,y=0x=1,y=0x=1,y=0。
当b!=0b!=0b!=0时，递归求解方程bx+(a%b)y=gcd(b,a%b)bx+(a\%b)y=gcd(b,a\%b)bx+(a%b)y=gcd(b,a%b)，直到b=0b=0b=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.青蛙的约会（怎么变蓝题了，我刚做的时候不是绿题吗QAQQAQQAQ）

---

大步小步算法（BSGS）

• 应用：求解高次同余方程
  ax≡b(modc)a^x \equiv b(mod\ c)ax≡b(mod c)
  其中gcd(a,c)=1gcd(a,c)=1gcd(a,c)=1。
• 由于我也证不来为什么算法是对的，所以我直接粘贴代码了qwqqwqqwq：
  看这篇大佬的博客吧：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】随机数生成器；
  ai−1≡xi+ba+1x1+a−1a^{i-1} \equiv \frac{x_i+\frac{b}{a+1}}{x_1+\frac{}{a-1}}ai−1≡x1​+a−1​xi​+a+1b​​
  其中xi=tx_i=txi​=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;
}

---

矩阵相关

• 矩阵
  这个就叫矩阵：
  [a11,a12,...,a1na21,a22,...,a2n...,...,...,...an1,an2,...,ann]\left[\begin{aligned} a_{11},a_{12},...,a_{1n}\\ a_{21},a_{22},...,a_{2n}\\ ...,...,...,...\\ a_{n1},a_{n2},...,a_{nn} \end{aligned}\right]⎣⎢⎢⎢⎢⎡​a11​,a12​,...,a1n​a21​,a22​,...,a2n​...,...,...,...an1​,an2​,...,ann​​⎦⎥⎥⎥⎥⎤​
• 矩阵运算

1.矩阵加减：
条件：两个矩阵长宽一样，例如
[142200]+[005750]=[1+04+02+52+70+50+0]=[147950]\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][1 4 22 0 0​]+[0 0 57 5 0​]=[1+0 4+0 2+52+7 0+5 0+0​]=[1 4 79 5 0​]
减法类似；
2.矩阵乘法：
条件：第一个矩阵的行数和第二个矩阵的列数相等。
[102−131]∗[312110]=[(1∗3+0∗2+2∗1)(1∗1+0∗1+2∗0)(−1∗3+3∗2+1∗1)(−1∗1+3∗1+1∗0)]=[5142]\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][1 0 2−1 3 1​]∗⎣⎢⎡​3 12 11 0​⎦⎥⎤​=[(1∗3+0∗2+2∗1) (1∗1+0∗1+2∗0)(−1∗3+3∗2+1∗1) (−1∗1+3∗1+1∗0)​]=[5 14 2​]
3.矩阵的幂:
[a11,a12,...,a1na21,a22,...,a2n...,...,...,...an1,an2,...,ann]2=[a11,a12,...,a1na21,a22,...,a2n...,...,...,...an1,an2,...,ann]∗[a11,a12,...,a1na21,a22,...,a2n...,...,...,...an1,an2,...,ann]\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]⎣⎢⎢⎢⎢⎡​a11​,a12​,...,a1n​a21​,a22​,...,a2n​...,...,...,...an1​,an2​,...,ann​​⎦⎥⎥⎥⎥⎤​2=⎣⎢⎢⎢⎢⎡​a11​,a12​,...,a1n​a21​,a22​,...,a2n​...,...,...,...an1​,an2​,...,ann​​⎦⎥⎥⎥⎥⎤​∗⎣⎢⎢⎢⎢⎡​a11​,a12​,...,a1n​a21​,a22​,...,a2n​...,...,...,...an1​,an2​,...,ann​​⎦⎥⎥⎥⎥⎤​

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

#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]=1a[1]=a[2]=a[3]=1a[1]=a[2]=a[3]=1
  a[x]=a[x−3]+a[x−1](x>3)a[x]=a[x-3]+a[x-1] (x>3)a[x]=a[x−3]+a[x−1](x>3)
  求aaa数列的第nnn项对1000000007（109+7）（10^9+7）（109+7）取余的值。
  nnn很大，如果递推求会超时。
  但是如果我们把这个模型放在矩阵中，求第二项就会变成求这个东西：
  [100010101]2\left[\begin{aligned}1\ 0\ 0\\0\ 1\ 0\\1\ 0\ 1 \end{aligned}\right]^2⎣⎢⎡​1 0 00 1 01 0 1​⎦⎥⎤​2
  求第nnn项：
  [100010101]n\left[\begin{aligned}1\ 0\ 0\\0\ 1\ 0\\1\ 0\ 1 \end{aligned}\right]^n⎣⎢⎡​1 0 00 1 01 0 1​⎦⎥⎤​n
  输出(1,1)(1,1)(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~xxx-1中与xxx互质的数的个数。
• 特点：欧拉函数是积性函数。
• 代码：

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;}//质数的phifor(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;
}

---

欧拉定理及费马小定理

• 欧拉定理：若aaa与mmm互质，则aφ(m)≡1(modm)a^{\varphi(m)}\equiv 1(mod\ m)aφ(m)≡1(mod m);
• 拓展欧拉定理：若r>φ(m)r>\varphi(m)r>φ(m)，则ar≡armodφ(m)+φ(m)(modm)a^r\equiv a^{r\ mod\ \varphi(m)+\varphi(m)}(mod\ m)ar≡ar mod φ(m)+φ(m)(mod m);
• 费马小定理：ap−1≡1(modq)a^{p-1}\equiv 1(mod\ q)ap−1≡1(mod q)，其中qqq为质数。（费马小定理是欧拉定理的特殊形式）
• 代码
  欧拉定理：

#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;
}

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

---

中国剩余定理

• 作用：用于求解同余方程组：
  {x≡b1(modm1)x≡b2(modm2)................x≡bn(modmn)\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.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​x≡b1​(mod m1​)x≡b2​(mod m2​)................x≡bn​(mod mn​)​
  其中m1,m2,...,mnm_1,m_2,...,m_nm1​,m2​,...,mn​互质。
• 代码

#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);elseprintf("No");return 0;
}

---

卢卡斯定理

• 作用
  求Cnm%pC^m_n\%pCnm​%p(ppp为质数)的值。
• 定理：
  1.Cnm=Cn/pm/p∗Cn%pm%pC^m_n=C^{m/p}_{n/p}*C^{m\%p}_{n\%p}Cnm​=Cn/pm/p​∗Cn%pm%p​
  2.将n,mn,mn,m写成ppp进制：
  Cnm=Cnk−1mk−1∗Cnk−2mk−2∗...∗Cn0m0(modp)C^m_n=C^{m_{k-1}}_{n_{k-1}}*C^{m_{k-2}}_{n_{k-2}}*...*C^{m_0}_{n_0}(mod\ p)Cnm​=Cnk−1​mk−1​​∗Cnk−2​mk−2​​∗...∗Cn0​m0​​(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;
}

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

---

拓展卢卡斯定理

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

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

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

---

莫比乌斯函数

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

#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;}//质数的phifor(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.∑d∣nμ(d)=[n=1]\sum_{d|n}\mu(d)=[n=1]d∣n∑​μ(d)=[n=1]
  也就是说，当n=1n=1n=1时，上式的值为1，否则为0.
  2.∑d∣nφ(d)=n\sum_{d|n}\varphi(d)=nd∣n∑​φ(d)=n
• 结论
  给定数论函数fff和ggg：
  若fff满足f(n)=∑d∣ng(d)f(n)=\sum_{d|n}g(d)f(n)=∑d∣n​g(d)，则ggg满足$ g(n)=\sum_{d|n}\mu(d)*f(\frac{n}{d})$
• 证明
  ∑d∣nμ(d)∗f(nd)=∑d∣nμ(d)∑d′∣ndg(d′)=∑d′∣ng(d′)∑d∣nd′μ(d)=g(n)\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)∑d∣n​μ(d)∗f(dn​)=∑d∣n​μ(d)∑d′∣dn​​g(d′)=∑d′∣n​g(d′)∑d∣d′n​​μ(d)=g(n)
• 应用
  将很复杂的式子变形成为很简单的式子。
• 例题
  1.【SDOI2015】约数个数和；
  2.【POI2007】ZAP-Queries；
  双倍经验：
  3.GCD；
  4.YY的GCD；
• 总结
  这种题通常是先把复杂的式子化为简单的式子，再用数论分块来做。
  这里以ZAP为例：
  目标公式：∑i=1m∑j=1n[gcd(i,j)=k]\sum^m_{i=1}\sum^n_{j=1}[gcd(i,j)=k]i=1∑m​j=1∑n​[gcd(i,j)=k]
  转化：
  ∑i=1m∑j=1n[gcd(i,j)=k]\sum^m_{i=1}\sum^n_{j=1}[gcd(i,j)=k]∑i=1m​∑j=1n​[gcd(i,j)=k]
  =∑i=1m∑j=1n∑d∣gcd(i,k)/kμ(d)=\sum^m_{i=1}\sum^n_{j=1}\sum_{d|gcd(i,k)/k}\mu(d)=∑i=1m​∑j=1n​∑d∣gcd(i,k)/k​μ(d)
  =∑dμ(d)∑kd∣im∑kd∣jn1=\sum_d\mu(d)\sum^m_{kd|i}\sum^n_{kd|j}1=∑d​μ(d)∑kd∣im​∑kd∣jn​1
  =∑dμ(d)⌊mkd⌋⌊nkd⌋=\sum_d\mu(d)\lfloor \frac{m}{kd} \rfloor \lfloor \frac{n}{kd} \rfloor=∑d​μ(d)⌊kdm​⌋⌊kdn​⌋;
  预处理⌊mkd⌋\lfloor \frac{m}{kd} \rfloor⌊kdm​⌋和⌊nkd⌋\lfloor \frac{n}{kd} \rfloor⌊kdn​⌋
• 代码：

#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;
}

---

杜教筛

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

即：
S(n)=∑i=1nh(i)−∑i=2ng(i)S(⌊ni⌋)g(1)S(n)=\frac{\sum^n_{i=1}h(i)-\sum^n_{i=2}g(i)S(\lfloor \frac{n}{i} \rfloor)}{g(1)}S(n)=g(1)∑i=1n​h(i)−∑i=2n​g(i)S(⌊in​⌋)​
用数论分块即可求出正解。
总复杂度:O(n34)O(n^{\frac{3}{4}})O(n43​)

• 模板
  先预处理出小范围的答案，大范围的用map存。
  求∑i=1nμ(i)\sum^n_{i=1} \mu(i)∑i=1n​μ(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;}//质数的phifor(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+bia+bia+bi的数，其中i2=−1i^2=-1i2=−1。
  2.四则运算：
  a.(a+bi)+(c+di)=(a+c)+(b+d)i(a+bi)+(c+di)=(a+c)+(b+d)i(a+bi)+(c+di)=(a+c)+(b+d)i;
  b.(a+bi)−(c+di)=(a−c)+(b−d)i(a+bi)-(c+di)=(a-c)+(b-d)i(a+bi)−(c+di)=(a−c)+(b−d)i;
  c.(a+bi)∗(c+di)=(ac−bd)+(bc+ad)i(a+bi)*(c+di)=(ac-bd)+(bc+ad)i(a+bi)∗(c+di)=(ac−bd)+(bc+ad)i;
  d.除法被吞了qwqqwqqwq;

• 单位复数根：
  1.欧拉公式:eix=cos⁡x+isin⁡xe^{ix}=\cos x+i\sin xeix=cosx+isinx；
  2.设ωnk\omega^k_nωnk​表示nnn的kkk次复根，即把一个平面均分成nnn份，取逆时针前kkk份。当n=8n=8n=8时，ωnk\omega^k_nωnk​如下：
  
  3.根据欧拉公式：ωnk=cos⁡(2πkn)+isin⁡(2πkn)\omega^k_n=\cos (\frac{2\pi k}{n})+i\sin (\frac{2\pi k}{n})ωnk​=cos(n2πk​)+isin(n2πk​)
  4.性质：
  a.ωnn=ωn0=1\omega^n_n=\omega^0_n=1ωnn​=ωn0​=1;
  b.消去引理：ωdndk=ωnk\omega^{dk}_{dn}=\omega^{k}_{n}ωdndk​=ωnk​;
  c.折半引理：(ωnk)2=ωn2k(\omega^k_n)^2=\omega^k_{\frac{n}{2}}(ωnk​)2=ω2n​k​;
  d.ωn1k1∗ωn2k2=ωn1+n2k1+k2\omega^{k_1}_{n_1}*\omega^{k_2}_{n_2}=\omega^{k_1+k_2}_{n_1+n_2}ωn1​k1​​∗ωn2​k2​​=ωn1​+n2​k1​+k2​​;

• 快速傅里叶变换
  终于进入正题了，写的我好辛苦啊！！！
  1.应用：用于求多项式的乘积，甚至可以求A∗BProblemA*B\ ProblemA∗B Problem(雾)。
  2.算法：Cooley-Tukey（分治）
  一般求多项式的乘积用的是系数表示法A=∑i=0nai∗xiA=\sum^n_{i=0}a_i*x^iA=∑i=0n​ai​∗xi，这样求多项式的乘积的时间复杂度就是O(n2)O(n^2)O(n2)，有点慢。
  而FFT选用的是点值表示法(x1,y1),...,(xn,yn)(x_1,y_1),...,(x_n,y_n)(x1​,y1​),...,(xn​,yn​)，多项式的乘积就是(x1,ya1∗yb1),...,(xn,yan∗ybn)(x_1,y_{a_1}*y_{b_1}),...,(x_n,y_{a_n}*y_{b_n})(x1​,ya1​​∗yb1​​),...,(xn​,yan​​∗ybn​​)，时间复杂度是O(n)O(n)O(n)，好快啊！！！
  接下来就是选择点值表达式中的点了。上面我们提到了单位复数根，它有很多优美的性质，我们就选它作为点横坐标啦(雾)。
  那么:A(ωnk)=∑i=0n−1ai∗ωnkiA(\omega^k_n)=\sum^{n-1}_{i=0}a_i*\omega^{k_i}_nA(ωnk​)=i=0∑n−1​ai​∗ωnki​​
  把A(x)A(x)A(x)的系数按奇偶分开，即：
  A0(x)=a0+a2x+a4x2+...A_0(x)=a_0+a_2x+a_4x^2+...A0​(x)=a0​+a2​x+a4​x2+...
  A1(x)=a1+a3x+a5x2+...A_1(x)=a_1+a_3x+a_5x^2+...A1​(x)=a1​+a3​x+a5​x2+...
  A(x)=A0(x2)+xA1(x2)A(x)=A_0(x^2)+xA_1(x^2)A(x)=A0​(x2)+xA1​(x2)
  有兴趣可以自己验证一下。
  那么：A(ωnk)=∑i=0n2−1a2iωn2ki+ωnk∑i=0n2−1a2i+1ωn2kiA(\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}A(ωnk​)=i=0∑2n​−1​a2i​ωn2ki​+ωnk​i=0∑2n​−1​a2i+1​ωn2ki​
  当k&lt;n2k&lt;\frac{n}{2}k<2n​时：
  A(ωnk)=∑i=0n2−1a2iωn2ki+ωnk∑i=0n2−1a2i+1ωn2kiA(\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(ωnk​)=i=0∑2n​−1​a2i​ω2n​ki​+ωnk​i=0∑2n​−1​a2i+1​ω2n​ki​
  A(ωnk+n2)=∑i=0n2−1a2iωn2ki+ωnk∑i=0n2−1a2i+1ωn2kiA(\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}}A(ωnk+2n​​)=i=0∑2n​−1​a2i​ω2n​ki​+ωnk​i=0∑2n​−1​a2i+1​ω2n​ki​
  每次代入规模减半，复杂度为O(nlog⁡n)O(n\log n)O(nlogn)

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

• 实现
  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)=∑i=1n−1a[i]+b[n−i]c(n)=\sum^{n-1}_{i=1}a[i]+b[n-i]c(n)=∑i=1n−1​a[i]+b[n−i]的式子适合用fftfftfft或nttnttntt解决。

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

---

快速数论变换（NTT）

• 原根
  若ggg是ppp的一个原根，则g0,g1,...,gp−1(modp)g^0,g^1,...,g^{p-1}(mod\ p)g0,g1,...,gp−1(mod p)是p−1p-1p−1个非0数。

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

• 代码
  改天再贴。

---

*斯特林数相关算法

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

算法笔记——数学相关相关推荐

1. 算法笔记-链相关、链的基础、单链双链环链、链的各种功能实现、链的算法题、面试题以及算法优化方法（多）、C#

   1. 链定义及其基础 单链表是一种链式存取的数据结构,用一组地址任意的存储单元存放线性表中的数据元素.这组存储单元既可以是连续的,也可以是不连续的. 链表定义: 链表是一种线性表数据结构: 从底层存储 ...

2. 算法笔记-堆相关、堆的定义、大小根堆、算法程序实现、堆的算法题、C#写法

   内容概述 1,堆结构就是用数组实现的完全二叉树结构 2,完全二叉树中如果每棵子树的最大值都在顶部就是大根堆 3,完全二叉树中如果每棵子树的最小值都在顶部就是小根堆 4,堆结构的heaplnsert与h ...

3. ACM模板 | 学习笔记 数学相关

   正在整理中qwq 欧拉函数 求欧拉函数的个数(phi(n)) phi(n)=n∗(1−1p1)∗...∗(1−1pk)phi(n) = n * (1 - \frac{1}{p_1})*...*(1-\ ...

4. 【 反向传播算法 Back-Propagation 数学推导以及源码详解 深度学习 Pytorch笔记 B站刘二大人（3/10）】

   反向传播算法 Back-Propagation 数学推导以及源码详解 深度学习 Pytorch笔记 B站刘二大人(3/10) 数学推导 BP算法 BP神经网络可以说机器学习的最基础网络.对于普通的简单 ...

5. 数学建模算法笔记（2）——主成分分析

   数学建模算法笔记(2)–主成分分析 目的:主成分分析的主要目的是希望用较少的变量去解释原来资料中的大部分变异,将我 们手中许多相关性很高的变量转化成彼此相互独立或不相关的变量,实际上是一种降维方法. ...

6. 《算法笔记》中文版 - 包括数组，链表，树，图，递归，DP，有序表等相关数据结构与算法的讲解及代码实现...

   来源:专知本文为资源,建议阅读5分钟本文为你分享<算法笔记>中文版. https://github.com/Dairongpeng/algorithm-note 目录概览 第一节 复杂度. ...

7. FE之DR之线性降维：PCA/白化、LDA算法的数学知识(协方差矩阵)、相关论文、算法骤、代码实现、案例应用等相关配图之详细攻略

   FE之DR之线性降维:PCA/白化.LDA算法的数学知识(协方差矩阵).相关论文.算法骤.代码实现.案例应用等相关配图之详细攻略 目录 PCA 1.PCA的数学知识 1.协方差矩阵计算 2.PCA算法 ...

8. 不明觉厉！Gitee大神们的算法/数学相关开源项目推荐

   现在的大厂面试,算法似乎已经成为了必考项目.当大家的业务水平相近,谁的数学与算法基础更好,谁可能就会获得更好的机会.Gitee 上也有一些数学算法的大牛,今天就为大家分享他们的开源项目,希望能给正在学 ...

9. Python学习笔记17：标准库之数学相关（math包，random包）

   前面几节看得真心累.如今先来点简单easy理解的内容. 一 math包 math包主要处理数学相关的运算. 常数 math.e   # 自然常数e math.pi  # 圆周率pi 运算函数 math ...

最新文章

1. html5--5-15 绘制阴影
2. 深入理解Java虚拟机（一）：Java内存模型
3. html 制作条形图,Highcharts 基本条形图
4. 深入浅出面向对象和原型【番外篇——重新认识new】
5. HarmonyOS应用开发——使用HUAWEI DevEco Studio创建第一个程序 HELLO WORLD！
6. java 微信支付实现
7. 前端学习（778）：随机数方法
8. 基于OOS批量修改资源标签值
9. 认识死锁之死锁的基本概念
10. 关于npm邮箱验证问题
11. 8个很棒的 jQuery 倒计时插件和教程
12. 得物App获得2020“年度新经济企业”奖项
13. 力扣题目——997. 找到小镇的法官
14. 315页 A Tutorial on Graph Neural Networks for NLP
15. 最新微信小程序抓包方法
16. java实现从浏览器读取Csv文件解析成 ListMap
17. HTML中的window对象和document对象详解
18. 仿真Windows_XP画图板的java实现
19. 产品全类目下找不到关键词，只有到特定类目才能找到，修改PRODUCT TYPE
20. 微信小程序项目实战：电影购票系统-李宁-专题视频课程

热门文章

1. Lora1278驱动V4.4.2讲解一：驱动移植
2. 开始面试，我该准备什么
3. 深圳可能是理解中国数字城市建设的窗口
4. NISP一级模拟题（一、二）
5. 命令查看spf_什么是SPF、邮箱域名SPF记录查询方法
6. 算法设计-利用栈判别表达式中的括弧是否配对
7. 完美国际真数苹果_苹果手机:Checkm8漏洞永久性破解A5-A11设备 全线旧设备实现完美越狱...
8. 使用Graham扫描法获取一个平面点集的凸包
9. win10 更新后摄像头问题
10. 怎么把字母缩小当符号_iPhone自带输入法怎么用 iPhone自带输入法小技巧【详解】...