神O的数论全家桶

§ (0.1) 前言

已经听 ~~神橡树~~ 欧稳欧讲了两天的数论了，最后一道题总是只能拿部分分，估计是掌握得不太熟练，写篇博客复习复习。

~~另外，某deco天天闹着自己要爆零，结果AK了两天的比赛，假人。下次不帮他点外卖~~

“这些都是水题！” ——deco这么说道，留下一脸无奈的tqr

~~那就删掉吧~~

§ (1.1) 扩展欧几里得算法

求 (ax+by=gcd(a,b)) 的一组解，可以用扩展欧几里得算法

设 ((x',y')) 是 (bx'+(a mod b)y'=gcd(a, b)) 的一组解

注意到

[a mod b=a-bleftlfloordfrac{a}{b} ight floor ]

对上式变形得到

[egin{cases}x=y'\y=x'-leftlfloordfrac{a}{b} ight floorend{cases} ]

时间复杂度同欧几里得算法，为 (Theta(log_{phi}a))

令 (d=gcd(a,b))

如果要求 (ax+by=c) 的解，由裴蜀定理 (d|c)，将扩展欧几里得算法的解乘 (dfrac{c}{d}) 即可

上面求出的只是 (ax+by=c) 的一组特解 ((x_0,y_0))

要求通解，可以设 (x=x_0+s，y=y_0-t)，带入原方程，得

[as=bt ]

即

[s=dfrac{b}{a}t=dfrac{b/d}{a/d}t ]

也就是说通解为（其中 (kinmathbb{Z})）

[egin{cases}x=x_0+frac{b}{d}k\y=y_0-frac{a}{d}kend{cases} ]

§ (1.2) 同余及其性质

若 (aequiv bpmod{m})，则

(apm cequiv bpm cpmod{m})
(acequiv bcpmod{m})
若 (c|a，c|b，则dfrac{a}{c}equivdfrac{b}{c} pmod{m/gcd(m,c) })

形如 (axequiv bpmod{m}) 的方程称为同余方程

根据同余的定义，同余方程等价于 (ax+mt=b (t in mathbb{Z}))，可以用扩展欧几里得求解

同余方程有解的条件是 (gcd(a,m) | b)

§ (1.3) 乘法逆元及求法

在题目需要取模的情况下，加，减，乘，都可以直接运算后取模，但除法不行

举个例子，(ans=ans \% mod ÷2)，直接除的话对答案会有影响，因此我们需要使用逆元（记作：(inv(n))）

(ans=ans \% mod ×inv(2)) 就可以了！

可以在 (Theta(n)) 内求出 (1) 到 (n) 的所有数模素数 (p) 的逆元，要求 (n<p)

首先 (1^{-1}=1)

对于一个数 (x)（(x>1)），设 (p=ax+b) ，

则 (a=leftlfloordfrac{p}{x} ight floor)，(b=p mod x)

有

[ax+bequiv 0pmod{p} ]

两边同时乘以 (b^{-1}x^{-1})，得：

[ab^{-1}+x^{-1}equiv0pmod{p} ]

即

[x^{-1}equiv-ab^{-1}equivleftlfloordfrac{p}{x} ight floor(p mod x)^{-1}pmod{p} ]

因为(p mod x<x)，所以其逆元已经计算过了，可以按 (x) 从小到大的顺序递推求解

以下是代码

inv[1]=1;
for(register int i=2;i<=n;++i)
{
	inv[i]=mod-mod/i*inv[mod%i]%mod;
	printf("%lld
",inv[i]);
}

§ (1.4) 中国剩余定理

设 (m_1,m_2,m_3,...,m_n) 两两互质，则同余方程组

[egin{cases} xequiv a_1pmod{m_1}\ xequiv a_2pmod{m_2}\ ......\ xequiv a_npmod{m_n}\ end{cases}]

在膜 (M=prod{m_i}) 下有唯一解

[xequiv sum{a_it_iM_i} ]

其中(M_i=M/m_i)，(t_i) 为 (M_i) 在膜 (m_i) 下的逆元

§ (1.5) 约数及约数相关性质

将 (n) 质因数分解

[n=prod{p_i^{r_i}} ]

约数个数

[sigma_0(n)=prod{(1+r_i)} ]

约数和

[sigma_1(n)=prodsumlimits_{k=0}^{r_i}p_i^{kt} ]

约数函数

[sigma_t(n)=sumlimits_{d|n}d^t=prodsumlimits_{k=0}^{r_i}p_i^{kt} ]

不大于 (sqrt{n}) 的约数不超过 (sqrt{n}) 个，大于 (sqrt{n}) 的约数 (d) 唯一对应一个小于 (sqrt{n}) 的约数 (dfrac{n}{d}) ，也不会超过 (sqrt{n}) 个，所以约数个数的上界是 (O(sqrt{n}))

随机数据下，约数个数的期望是 (O(ln_{ n}))

§ (1.6) 欧拉函数与(扩展)欧拉定理

上面我们已经知道，欧拉函数 (phi(n)) 表示不超过 (n) 与 (n) 互质的数的个数

欧拉函数可以用莫比乌斯函数容斥，枚举公约数 (d)

[phi(n)=sumlimits_{d|n}mu(d)dfrac{n}{d} ]

将 (n) 质因数分解

[n=prod{p_i^{r^i}} ]

代入莫比乌斯函数定义式，得

[phi(n)=nprod(1-dfrac{1}{p_i})=prod{p_i^{r_i-1}}(p_i-1) ]

使用线性筛 (Theta(n)) 地求出 (1) 到 (n) 所有数的欧拉函数值，由表达式可以发现，欧拉函数也有积性

对于质数 (p)，(phi(p)=p-1)

对于 (px)，如果 (p mid x)，(phi(p)phi(x)=(p-1)phi(x))

对于 (px)，如果 (p|x)，(phi(px)=pphi(x))

费马小定理

对于素数 (p) 和任意正整数 (ain[1,p))，有

[a^{p-1}equiv 1pmod{p} ]

事实上，费马小定理是欧拉定理的特殊情况

欧拉定理

给定正整数 (a,m)，若 (gcd(n,m)=1)，则

[a^{phi(m)}equiv 1pmod{m} ]

扩展欧拉定理

给定正整数 (a,m)，(r>phi(m))，则

[a^requiv a^{r mod phi(m)+phi(m)}pmod{m} ]

欧拉定理常用来求逆元

[a^{-1}equiv a^{phi(m)-1}pmod{m} ]

特别的，当 (m) 是素数 (p) 时，

[a^{-1}equiv a^{p-2}pmod{p} ]

使用快速幂实现，时间复杂度为 (Theta(lg{m}))

§ (1.7) BSGS算法及扩展

给定正整数 (a,b,m)，求最小的非负整数 (x)，满足

[a^xequiv bpmod{m} ]

显然，由欧拉定理，(x<phi(m))

但如果 (gcd(a,m)=1)，可以用BSGS算法解决

先特判 (x=0)，否则，选定一个参数 (s)，设 (x=is-j)，其中 (0leq j<s)

变形可得

[a^{is}equiv a^jbpmod{m} ]

从小到大枚举 (i)，将 (a^{is} mod m) 存入哈希表，相同的保留最小的 (i)，复杂度为 (Theta(phi(m)/s))

然后从小到大枚举 (j)，在哈希表中查询 (a^jb)，这一步复杂度为 (Theta(s))

BSGS算法可以处理给定 (a,m,q) 次询问给出 (b) 的问题，复杂度为 (Theta(phi(m)/s+qs))，取 (s=sqrt{phi(m)/q}) 时最优，复杂度为 (Theta(sqrt{qphi(m)}))

扩展BSGS算法

如果 (a,m) 不互质，就需要进行一些扩展

如果 (b=0)，则 (m=1) 时有解，否则无解

如果 (b=1)，则 (x=0)

否则，设 (d=gcd(a,m))，如果 (d mid b) 无解，否则两边同除以 (d)，则

[(a/d)a^{x-1}equiv b/dpmod{m/d} ]

因为 (a/d，m/d) 互质

[a^{x-1}equiv(a/d)^{-1}(b/d)pmod{m/d} ]

多次执行上面的过程，直到 (a,m) 互质，然后使用BSGS求解

§ (1.8) 组合数取模

求

[C_m^n mod p ]

如果 (n,m) 较小，可以 (Theta(nm)) 递推（杨辉三角了解一下）

[C_m^n=C_{m-1}^n+C_{m-1}^{n-1} ]

如果 (n,m) 不太大，且 (p) 是质数，可以 (Theta(n)) 预处理阶乘及阶乘逆元

[C_m^n=dfrac{n!}{m!(n-m)!} ]

如果 (n,m) 很大，(p) 是质数且不太大，可以用卢卡斯定理，将 (n,m) 写成 (p) 进制

[n=sumlimits_{i=0}^{l-1}n_ip^i，m=sumlimits_{i=0}^{l-1}m_ip^i ]

[C_m^n=C_{m_{k-1}}^{n_{k-1}}*C_{m_{k-2}}^{n_{k-2}}*…*C_{m_0}^{n_0} mod p ]

需要预处理 (0) 到 (p-1) 的阶乘及阶乘逆元，复杂度为预处理 (Theta(p))，单次询问 (Theta(log_{p}n))

如果 (p) 不是质数该怎么办呢？

将 (p) 质因数分解

[p=prod{p_i^{r_i}} ]

对每个 (p^r) 分别取膜计算，然后用中国剩余定理合并

将 (n!,m!,(n-m)!) 表示为 (ap^i)，其中 (p mid a)，这样就可以用 (a) 的逆元，(p) 的指数相加减即可

设答案为 (f(n))，将 (n!) 中所有 (p) 的倍数提出来，这一部分是 (f(lfloor{n/p} floor)p^{lfloor{n/p} floor})

预处理 (0) 到 (p^r-1) 去除 (p) 的倍数的阶乘 (g(i))

则剩下的部分以 (p^r) 为循环节，有

[f(n)=g^{lfloor{n/p^r} floor}(p^r-1)g(n mod p^r)f(lfloor{n/p} floor)p^{lfloor{n/p} floor} ]

§ (2.1) 一些基本函数

数论函数

定义域为正整数的函数成为数论函数

积性函数

对于数论函数 (f)，若任意互质的 (p，q) 都有 (f(pq)=f(p)f(q))，则称 (f) 是积性函数

完全积性函数

对于数论函数 (f)，若任意 (p，q) 都有 (f(pq)=f(p)f(q))，则称 (f) 是完全积性函数

常见积性函数

除数函数 (sigma_k(n))：(n)的所有约数的 (k) 次方和

欧拉函数 (phi(n))：不超过(n)与(n)互质的数的个数

莫比乌斯函数 (mu(n)=egin{cases}0&exists p,p^2|n\(-1)^r&n=p_1p_2…end{cases})

（常用于约数相关的容斥）

常见完全积性函数

1函数 (1(n)=1)

原函数 (varepsilon(n)=[n=1])

幂函数 (I^k(n)=n^k)

§ (2.2) 狄利克雷卷积¤

更新中

§ (2.3) 莫比乌斯反演¤

更新中

§ (2.4) 杜教筛¤

更新中

§ (2.5) 洲阁筛¤

更新中

§ (3.1) FFT快速傅里叶变换

多项式的表示法

系数表示法：

设 (A(x)) 表示一个 (n-1) 次多项式

则 (A(x)=sum_{i=0}^{n}a_i*x^i)

看着很复杂，但其实就是一个标准多项式

例如 (A(3)=x^2+3x+2)

使用这种方法，计算多项式乘法的复杂度为 (Theta(n^2))
点值表示法

如果将 (n) 个互不相同的 (x) 代入多项式，会得到 (n) 个不同取值的 (y)

类似于两点确定一条直线（即一次多项式），这个 (n-1) 次多项式被这 (n) 个点 ((x_1,y_1),(x_2,y_2),...,(x_n,y_n)) 唯一确定

其中 (y_i=sum_{j=0}^{n-1}a_j*x^i)

例如上面的 (A(3)=x^2+3x+2) ，就可以表示为 ((0,2),(1,5),(2,12))

这种方法计算多项式乘法的时间复杂度仍然为 (Theta(n^2))

可以看到，两种表示方法的时间复杂度相同（且很大），因此我们需要想办法对其进行优化

对于系数表示法，由于每次项的系数是固定的，因此优化困难

对于点值表示法，由于我们可以任意地选取不同的点，因此我们可以想办法尽可能地选取一些特殊的点，使其计算简便

复数的引入

如果您不明白向量是什么，不保证能理解此部分内容

定义

设 (a,b) 为实数， (i^2=-1)，则形如 (a+bi) 的数交复数，其中 (i) 是虚数单位，复数域是目前已知最大的域

在复平面中， (x) 轴代表实数， (y) 轴代表虚数，从原点 ((0,0)) 到 ((a,b)) 的向量表示复数 (a+bi)
- 模长
  从原点 ((0,0)) 到点 ((a,b)) 的距离，即 (sqrt{a^2+b^2})
- 幅角
  以逆时针为正方向，从 (x) 轴正半轴到已知向量的转角的有向角叫做幅角
运算法则
- 加法
  
  在复平面中，复数加法与向量加法相同，均满足平行四边形性质（(overrightarrow{AB}+overrightarrow{AD}=overrightarrow{AC})）
- 乘法
  
  几何定义：复数相乘，模长相乘，幅角相加
  
  代数定义：
  
  [(a+bi)*(c+di) ]
  [=ac+adi+bci+bdi^2 ]
  [=ac+adi+bci-bd ]
  [=(ac-bd)+(bc+ad) ]
单位根

下文中默认 (n) 为 (2) 的正整数次幂

在复平面上，以原点为圆心， (1) 为半径作圆，所得的圆叫单位圆。以圆心为起点，圆的 (n) 等分点为终点，作 (n) 个向量，设幅角为正且最小的向量对应的复数为 (omega_n)，称为 (n) 次单位根

根据复数乘法的运算法则，其余 (n-1) 个复数为 (omega_n^2,omega_n^3,...,omega_n^n)

其中 (omega_n^0=omega_n^n=1)

由欧拉公式，

[omega_n^k=cos{k*dfrac{2pi}{n}}+isin{k}*dfrac{2pi}{n} ]

单位根为幅角的 (dfrac{1}{n})

在单位根中，若 (z^n=1)，我们把 (z) 称为 (n) 次单位根

单位根的性质

[omega_{n}^{frac{n}{2}}=cos{dfrac{n}{2}}*dfrac{2pi}{n}+isin{dfrac{n}{2}*dfrac{2pi}{n}} ]

[=cos{pi}+isin{pi} ]

[=-1 ]

快速傅里叶变换

因为一个 (n) 次多项式可以由 (n) 个不同的点唯一确定

于是我们可以把单位根的 (0) 到 (n-1) 次幂代入，确定出多项式，但复杂度仍为 (Theta{(n^2)})

我们设多项式 (A(x)) 的系数为 ((a_0,a_1,...,a_{n-1}))

那么

[A(x)=a_0+a_1*x+a_2*x^2+a_3*x^3+...+a_{n-2}*x^{n-2}+a_{n-1}*x^{n-1} ]

将下标按奇偶分类，设

[A_1(x)=a_0+a_2*x+a_4*x^2+...+a_{n-2}*x^{frac{n}{2}-1} ]

[A_2(x)=a_1*x+a_3*x+a_5*x^2+...+a_{n-1}*x^{frac{n}{2}-1} ]

那么

[A(x)=A_1(x^2)+xA_2(x^2) ]

将 (omega_n^k<frac{n}{2}) 代入得

[A(omega_n^k)=A_1(omega_n^{2k})+omega_n^kA_2(omega_{frac{n}{2}}^k) ]

[=A_1(omega_{frac{n}{2}}^k)+omega_n^kA_2(omega^k_{frac{n}{2}}) ]

同理，将 (omega_n^{k+frac{n}{2}}) 代入得

[A(omega_n^{k+frac{n}{2}}) ]

[=A_1(omega_n^{2k+n})+omega_n^{k+frac{n}{2}}(omega_{n}^{2k+n}) ]

[=A_1(omega_n^{2k}*omega_n^n)-omega_n^kA_2(omega_n^{2k}*omega_n^n) ]

[=A_1(omega_n^{2k})-omega_n^kA_2(omega_n^{2k}) ]

这两个式子只有常数项有差异！

因此我们枚举第一个式子的时候，可以 (Theta(1)) 地得到第二个式子的值

又因为第一个式子的 (k) 在取遍 (left[0,frac{n}{2}-1 ight]) 的时候， (k+frac{n}{2}) 取遍了 (left[frac{n}{2},n-1 ight])

所以我们把问题的规模缩小了一半！

分治以后递归求解！

时间复杂度 (Theta(nlog{n}))

快速傅里叶逆变换

然而事情还没完

在在日常生活中我们并不常用点值表达式

因此我们还需要把点值表示法转换为系数表示法，这就是傅里叶逆变换

((y_0,y_1,y_2,...,y_{n-1})) 为 (a_0,a_1,a_2,...,a_{n-1}) 的傅里叶变换

设另一个向量 ((c_0,c_1,c_2,...,c_{n-1})) 满足

[c_k=sumlimits_{i=0}^{n-1}y_i(omega_n^{-k})^i ]

即多项式 (B(x)=y_0,y_1x,y_2x^2,...,y_{n-1}x^{n-1}) 在 (omega_n^0,omega_n^{-1},omega_n^{-2},...,omega_{n-1}^{-(n-1)}) 处的点值表示

我们推导一下：

((c_0,c_1,c_2,...,c_{n-1})) 满足

[c_k=sumlimits_{i=0}^{n-1}y_i(omega_n^{-k})^i ]

[=sumlimits_{i=0}^{n-1}(sumlimits_{j=0}^{n-1}a_j(omega_n^i)^j)(omega_n^{-k})^i ]

[=sumlimits_{i=0}^{n-1}(sumlimits_{j=0}^{n-1}a_j(omega_n^j)^i)(omega_n^{-k})^i ]

[=sumlimits_{i=0}^{n-1}sumlimits_{j=0}^{n-1}a_j(omega_n^j)^i(omega_n^{-k})^i ]

[=sumlimits_{i=0}^{n-1}sumlimits_{j=0}^{n-1}a_j(omega_n^{j-k})^i ]

[=sumlimits_{j=0}^{n-1}a_j(sumlimits_{i=0}^{n-1}(omega_n^{j-k})^i) ]

设 (S(x)=sum_{i=0}^{n-1}x^i)

将 (omega_n^k) 代入得

[S(omega_n^k)=1+(omega_n^k)+(omega_n^k)^2+...+(omega_n^k)^{n-1} ]

当 (k !=0) 时

等式两边同乘 (omega_n^k) 得

[omega_n^kS(omega_n^k)=omega_n^k+(omega_n^k)^2+(omega_n^k)^3+...+(omega_n^k)^n ]

两式相减得

[omega_n^kS(omega_n^k)-S(omega_n^k)=(omega_n^k)^n-1 ]

[S(omega_n^k)=dfrac{(omega_n^k)^k-1}{omega_n^k-1} ]

[S(omega_n^k)=dfrac{(omega_n^n)^k-1}{omega_n^k-1} ]

[S(omega_n^k)=dfrac{1-1}{omega_n^k-1} ]

显然这个式子的分母不为 (0)，但分子为 (0)

因此，当 (K !=0) 时，(S(omega_n^k)=0)

(k=0) 时，(S(omega_n^k)=0)

对于刚才的式子

[c_k=sumlimits_{j=0}^{n-1}a_j(sumlimits_{i=0}^{n-1}(omega_n^{j-k})^i) ]

当 (j !=k) 时，(c_k=0)

当 (j=k) 时，值为 (n)

即

[egin{cases}c_k=na_k\a_k=dfrac{c_k}{n}\end{cases} ]

就可以得到点值与系数之间的关系

总结及优化

FFT其实就是把系数表示法转换为点值表示法，再转换为系数表示法

然而，虽然上面的思路正确，如果你使用递归的方式实现FFT的话，常数会非常大，因此我们需要优化！

观察原序列和反转后的序列，我们发现所需要的序列就是原序列下标的二进制翻转！

因此我们对序列按照下标奇偶分类没有任何意义

我们可以 (Theta(n)) 地利用某种操作得到我们要求的序列，然后不断地向上合并

代码如下：

#include<bits/stdc++.h>
#define PI 3.1415926535897932384626
#define N 10000005
using namespace std;

int n,m,limit=1;
int l,r[N];

struct Complex
{
	double x,y;
	Complex(double xx=0,double yy=0){x=xx;y=yy;}
}a[N],b[N];

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

inline void FFT(Complex *A,int type)
{
	for(register int i=0;i<limit;++i)
		if(i<r[i]) swap(A[i],A[r[i]]);//求出要迭代的序列
	for(register int mid=1;mid<limit;mid<<=1)
	{
		Complex Wn (cos(PI/mid),type*sin(PI/mid));
		for(register int R=mid<<1,j=0;j<limit;j+=R)//R是区间右端点，j表示已经到哪个位置
		{
			Complex w(1,0);//幂
			for(register int k=0;k<mid;++k,w=w*Wn)//枚举左半部分
			{
				Complex x=A[j+k],y=w*A[j+mid+k];
				A[j+k]=x+y;
				A[j+mid+k]=x-y;
			}
		}
	}
	return;
}

int main()
{
	scanf("%d%d",&n,&m);
	for(register int i=0;i<=n;++i) scanf("%lf",&a[i].x);
	for(register int i=0;i<=m;++i) scanf("%lf",&b[i].x);
	while(limit<=n+m) limit<<=1,++l;
	for(register int i=0;i<limit;++i)
		r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	//i是i/2左移一位，反转后就应该右移一位
	FFT(a,1);FFT(b,1);
	for(register int i=0;i<=limit;++i) a[i]=a[i]*b[i];
	FFT(a,-1);
	for(register int i=0;i<=n+m;++i) printf("%d ",(int)(a[i].x/limit+0.1));
	return 0;
}

至此，FFT就完全结束了

§ (3.2) NTT快速数论变换¤

更新中

§ (3.3) FWT快速沃尔什变换¤

更新中

§ (3.4) 斯特林数¤

第一类斯特林数

更新中

第二类斯特林数

更新中

§ (3.5) 斯特林反演¤

更新中

§ (0.2) 说明

带“¤”的就是正在码字的部分……

内容实在太多了

因为要考(NOIP(的后身))，所以省选内容先不写了

§ (0.3) 最后的话

~~欧稳欧太强了！~~

~~欧稳欧 AK IOI 预定！~~

~~那就删掉吧~~