神O的数论全家桶

目录

• 神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)\) 的一组解

注意到

\[a\ mod\ b=a-b\left\lfloor\dfrac{a}{b}\right\rfloor \]

对上式变形得到

\[\begin{cases}x=y'\\y=x'-\left\lfloor\dfrac{a}{b}\right\rfloor\end{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 \]

也就是说通解为（其中 \(k\in\mathbb{Z}\)）

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

§ \(1.2\) 同余及其性质

若 \(a\equiv b\pmod{m}\)，则

• \(a\pm c\equiv b\pm c\pmod{m}\)
• \(ac\equiv bc\pmod{m}\)
• 若 \(c|a，c|b，则\dfrac{a}{c}\equiv\dfrac{b}{c} \pmod{m/gcd(m,c)\ }\)

形如 \(ax\equiv b\pmod{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=\left\lfloor\dfrac{p}{x}\right\rfloor\)，\(b=p\ mod\ x\)

有

\[ax+b\equiv 0\pmod{p} \]

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

\[ab^{-1}+x^{-1}\equiv0\pmod{p} \]

即

\[x^{-1}\equiv-ab^{-1}\equiv\left\lfloor\dfrac{p}{x}\right\rfloor(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\n",inv[i]);
}

§ \(1.4\) 中国剩余定理

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

\[\begin{cases} x\equiv a_1\pmod{m_1}\\ x\equiv a_2\pmod{m_2}\\ ......\\ x\equiv a_n\pmod{m_n}\\ \end{cases}\]

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

\[x\equiv \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)=\prod\sum\limits_{k=0}^{r_i}p_i^{kt} \]

约数函数

\[\sigma_t(n)=\sum\limits_{d|n}d^t=\prod\sum\limits_{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)=\sum\limits_{d|n}\mu(d)\dfrac{n}{d} \]

将 \(n\) 质因数分解

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

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

\[\phi(n)=n\prod(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\nmid x\)，\(\phi(p)\phi(x)=(p-1)\phi(x)\)

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

费马小定理

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

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

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

欧拉定理

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

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

扩展欧拉定理

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

\[a^r\equiv 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^x\equiv b\pmod{m} \]

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

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

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

变形可得

\[a^{is}\equiv a^jb\pmod{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{q\phi(m)})\)

扩展BSGS算法

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

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

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

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

\[(a/d)a^{x-1}\equiv b/d\pmod{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=\sum\limits_{i=0}^{l-1}n_ip^i，m=\sum\limits_{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\nmid a\)，这样就可以用 \(a\) 的逆元，\(p\) 的指数相加减即可

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

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

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

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

§ \(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)=\begin{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{2\pi}{n}}+i\sin{k}*\dfrac{2\pi}{n} \]

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

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

• 单位根的性质

\[\omega_{n}^{\frac{n}{2}}=\cos{\dfrac{n}{2}}*\dfrac{2\pi}{n}+i\sin{\dfrac{n}{2}*\dfrac{2\pi}{n}} \]

\[=\cos{\pi}+i\sin{\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\right]\) 的时候， \(k+\frac{n}{2}\) 取遍了 \(\left[\frac{n}{2},n-1\right]\)

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

分治以后递归求解！

时间复杂度 \(\Theta(n\log{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=\sum\limits_{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=\sum\limits_{i=0}^{n-1}y_i(\omega_n^{-k})^i \]

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

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

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

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

\[=\sum\limits_{j=0}^{n-1}a_j(\sum\limits_{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=\sum\limits_{j=0}^{n-1}a_j(\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i) \]

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

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

即

\[\begin{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 预定！

那就删掉吧