• 神O的数论全家桶


    神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 预定!

    那就删掉吧

  • 相关阅读:
    委托
    apply()和call()
    Sql小技巧
    plsql中文乱码
    Windows8中使用IE8等低版本浏览器
    React Native
    谷歌浏览器添加flash白名单
    jsonp原理详解
    垂直居中
    window.moveTo(),window.moveBy()不生效
  • 原文地址:https://www.cnblogs.com/tqr06/p/11272979.html
Copyright © 2020-2023  润新知