- 神O的数论全家桶
- § $0.1$ 前言
- § $1.1$ 扩展欧几里得算法
- § $1.2$ 同余及其性质
- § $1.3$ 乘法逆元及求法
- § $1.4$ 中国剩余定理
- § $1.5$ 约数及约数相关性质
- § $1.6$ 欧拉函数与(扩展)欧拉定理
- § $1.7$ BSGS算法及扩展
- § $1.8$ 组合数取模
- § $2.1$ 一些基本函数
- § $2.2$ 狄利克雷卷积¤
- § $2.3$ 莫比乌斯反演¤
- § $2.4$ 杜教筛¤
- § $2.5$ 洲阁筛¤
- § $3.1$ FFT快速傅里叶变换
- § $3.2$ NTT快速数论变换¤
- § $3.3$ FWT快速沃尔什变换¤
- § $3.4$ 斯特林数¤
- § $3.5$ 斯特林反演¤
- § $0.2$ 说明
- § $0.3$ 最后的话
神O的数论全家桶
§ (0.1) 前言
已经听 神橡树 欧稳欧 讲了两天的数论了，最后一道题总是只能拿部分分，估计是掌握得不太熟练，写篇博客复习复习。
另外，某deco天天闹着自己要爆零，结果AK了两天的比赛，假人。下次不帮他点外卖
“这些都是水题！” ——deco这么说道，留下一脸无奈的tqr
那就删掉吧
§ (1.1) 扩展欧几里得算法
求 (ax+by=gcd(a,b)) 的一组解，可以用扩展欧几里得算法
设 ((x',y')) 是 (bx'+(a mod b)y'=gcd(a, b)) 的一组解
注意到
对上式变形得到
时间复杂度同欧几里得算法，为 (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)，带入原方程，得
即
也就是说通解为（其中 (kinmathbb{Z})）
§ (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)
有
两边同时乘以 (b^{-1}x^{-1})，得：
即
因为(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) 两两互质，则同余方程组
在膜 (M=prod{m_i}) 下有唯一解
其中(M_i=M/m_i)，(t_i) 为 (M_i) 在膜 (m_i) 下的逆元
§ (1.5) 约数及约数相关性质
将 (n) 质因数分解
约数个数
约数和
约数函数
不大于 (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)
将 (n) 质因数分解
代入莫比乌斯函数定义式，得
使用线性筛 (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,m)，若 (gcd(n,m)=1)，则
扩展欧拉定理
给定正整数 (a,m)，(r>phi(m))，则
欧拉定理常用来求逆元
特别的，当 (m) 是素数 (p) 时，
使用快速幂实现，时间复杂度为 (Theta(lg{m}))
§ (1.7) BSGS算法及扩展
给定正整数 (a,b,m)，求最小的非负整数 (x)，满足
显然，由欧拉定理，(x<phi(m))
但如果 (gcd(a,m)=1)，可以用BSGS算法解决
先特判 (x=0)，否则，选定一个参数 (s)，设 (x=is-j)，其中 (0leq j<s)
变形可得
从小到大枚举 (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，m/d) 互质
多次执行上面的过程，直到 (a,m) 互质，然后使用BSGS求解
§ (1.8) 组合数取模
求
如果 (n,m) 较小，可以 (Theta(nm)) 递推（杨辉三角了解一下）
如果 (n,m) 不太大，且 (p) 是质数，可以 (Theta(n)) 预处理阶乘及阶乘逆元
如果 (n,m) 很大，(p) 是质数且不太大，可以用卢卡斯定理，将 (n,m) 写成 (p) 进制
需要预处理 (0) 到 (p-1) 的阶乘及阶乘逆元，复杂度为预处理 (Theta(p))，单次询问 (Theta(log_{p}n))
如果 (p) 不是质数该怎么办呢？
将 (p) 质因数分解
对每个 (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) 为循环节，有
§ (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)
由欧拉公式，
单位根为幅角的 (dfrac{1}{n})
在单位根中，若 (z^n=1)，我们把 (z) 称为 (n) 次单位根
- 单位根的性质
快速傅里叶变换
因为一个 (n) 次多项式可以由 (n) 个不同的点唯一确定
于是我们可以把单位根的 (0) 到 (n-1) 次幂代入，确定出多项式，但复杂度仍为 (Theta{(n^2)})
我们设多项式 (A(x)) 的系数为 ((a_0,a_1,...,a_{n-1}))
那么
将下标按奇偶分类，设
那么
将 (omega_n^k<frac{n}{2}) 代入得
同理，将 (omega_n^{k+frac{n}{2}}) 代入得
这两个式子只有常数项有差异！
因此我们枚举第一个式子的时候，可以 (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})) 满足
即多项式 (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})) 满足
设 (S(x)=sum_{i=0}^{n-1}x^i)
将 (omega_n^k) 代入得
当 (k !=0) 时
等式两边同乘 (omega_n^k) 得
两式相减得
显然这个式子的分母不为 (0)，但分子为 (0)
因此，当 (K !=0) 时，(S(omega_n^k)=0)
(k=0) 时，(S(omega_n^k)=0)
对于刚才的式子
当 (j !=k) 时，(c_k=0)
当 (j=k) 时，值为 (n)
即
就可以得到点值与系数之间的关系
总结及优化
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 预定！
那就删掉吧