- 神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 预定!
那就删掉吧