数学小结(持续更新)
ez_lcw
·
2020-08-25 21:47:26
·
个人记录
数学小结
前言:这篇文章专门整理了各种有关数学的东西(包括原理和例题),主要是怕自己忘了。
整理的不齐全,而且写的都是比较基础的东西。
UPD:把代码和例题全删了。
@[toc]
扩展欧几里得(exgcd)
问方程 ax+by=c 的 x 的最小非负整数解。
根据贝祖定理,若 \gcd(a,b)|c ,则方程有解,否则无解。
那么我们先可以求出 ax+by=\gcd(a,b) 的一组解,然后 x 、y 同时乘上 \dfrac{c}{\gcd(a,b)} ,就可以得到 ax+by=c 的解了。至于怎么求 x 的最小非负整数解,我们先要将原方程约分,即 a 、b 、c 都除以 \gcd(a,b) ,然后再根据不定方程,让 x 模上 b 就好了。
现在的问题就是如何求出 ax+by=\gcd(a,b) 的一组解。
考虑用类似求 \gcd 的辗转相除法求解。
假设我们已经求出了 a'x'+b'y'=\gcd(a',b') 的一组解,其中 a'=b ,b'=a\bmod b=a-\lfloor\dfrac{a}{b}\rfloor\times b ,现在考虑如何求 x 、y 的一组解。
显然,根据辗转相除法,\gcd(a,b)=\gcd(b,a\bmod b)=\gcd(a',b') 。
那么有:
\begin{aligned}
a'x'+b'y'&=\gcd(a',b')\\
bx'+\left(a-\lfloor\dfrac{a}{b}\rfloor\times b\right)y'&=\gcd(a,b)\\
ay'+b\left(x'-\lfloor\dfrac{a}{b}\rfloor\times y'\right)&=\gcd(a,b)
\end{aligned}
可以推得:
\begin{cases}
x=y'\\
y=x'-\lfloor\dfrac{a}{b}\rfloor\times y'
\end{cases}
这样就可以求出 x 、y 的一组解了。
考虑边界情况:当 b=0 时,有 ax=\gcd(a,b)=a ,得 x=1 ,y=0 。
然后根据一开始说的方法就能求出 x 的最小非负整数解了。
求逆元
若 a\times x\equiv 1 \pmod p ,则称 x 为 a 在模 p 意义下的逆元,记作 a^{-1} ,通常来算 \dfrac{1}{a}\bmod p 的值。(保证 \gcd(a,p)=1 )
若 p\in prime ,根据费马小定理,a^{p-1}\equiv 1\pmod p ,即 a\times a^{p-2}\equiv 1 \pmod p ,那么答案即为 a^{p-2}\bmod p 。
若 p\not\in prime ,那么 a\times x\equiv 1 \pmod p 可以看做是 ax+py=1 的一组解。
中国剩余定理(crt)
问以下方程的最小非负整数解:
\begin{cases}
x\equiv a_1 \pmod{p_{1}}\\
x\equiv a_2 \pmod{p_{2}}\\
\cdots\\
x\equiv a_k \pmod{p_{k}}
\end{cases}
其中 p_1,p_2,\cdots,p_k 为两两互质的数。
而这个定理貌似就是在构造一个特殊解。
设 M=\prod\limits_{i=1}^{k}p_i ,m_i=\dfrac{M}{p_i} ,t_i 为 m_i 在模 p_i 意义下的逆元,则有 m_i\times t_i\equiv 1 \pmod{p_i} 。
那么这个方程的一个特殊解就是:x_0=\sum\limits_{i=1}^{k}a_im_it_i 。
证明:
对于任意的 i ,我们考虑 x_0\equiv a_i \pmod{p_{i}} 是否成立。
将 x_0 拆解成两部分之和:x_0=\sum\limits_{j=1,j\neq i}^{k}a_jm_jt_j+a_im_it_i 。
显然,对于任意的 j\neq i ,a_jm_jt_j\bmod p_i 肯定为 0 ,因为 m_j 的定义就是除了 p_j 之外的 p_i 之积,得 m_j\bmod p_i 为 0 。
那么就可以得到 \sum\limits_{j=1,j\neq i}^{k}a_jm_jt_j\equiv 0 \pmod{p_i} 。
然后因为 m_it_i\equiv 1 \pmod{p_i} ,得 a_im_it_i\equiv a_i \pmod{p_i} 。
所以 x_0=\sum\limits_{j=1,j\neq i}^{k}a_jm_jt_j+a_im_it_i\equiv a_i \pmod{p_{i}} 。
那么 x_0 就是 x 的一个特殊解了。
证毕。
然后又因为 x 一定是 aM+b 的形式,所以 x_{\min}=x_0\bmod M 。
扩展中国剩余定理(excrt)
问以下方程的最小非负整数解:
\begin{cases}
x\equiv a_1 \pmod{p_{1}}\\
x\equiv a_2 \pmod{p_{2}}\\
\cdots\\
x\equiv a_k \pmod{p_{k}}
\end{cases}
发现这类问题和普通中国剩余定理的区别就是少了这个条件:“其中 p_1,p_2,\cdots,p_k 为两两互质的数”。
这就导致了 m_i 在模 p_i 意义下的逆元可能根本不存在,所以不能那么做。
假设前 i 个方程组的特殊解已经求出来为 x_0 。
设 M_i=\prod\limits_{j=1}^{i}p_j 。
那么前 i 个方程组的通解为 x_i=x_0+M_i\times t (t\in Z )。
我们现在就是要找到一个 t ,使得 x_0+M_i\times t\equiv a_i\pmod{p_{i}} ,即 M_i\times t\equiv a_i-x_0\pmod{p_{i}} 。
这个用 exgcd 来做就可以了。
然后这么一直下去就能得到最后的解了。
卢卡斯定理(Lucas)
问 \dbinom{n}{m}\bmod p 的值,其中 p\in prime 。
设 n=ap+b ,m=lp+r 。
那么 \dbinom{n}{m} 就相当于 (x+1)^n 的展开式中 x^m 的系数。
然后先证明一个东西:(x+1)^p\equiv x^p+1\pmod p 。
证明比较简单,把 (x+1)^p 展开,得 (x+1)^p\equiv \sum\limits_{i=0}^{p}\dbinom{p}{i}x^i\pmod p ,其中对于任意的 0<i<p ,肯定有 \dbinom{p}{i}=\dfrac{p\times (p-1)\times \cdots\times(p-i+1)}{1\times 2\times \cdots\times i} 是 p 的倍数,因为 p 为质数,1\sim i 均与 p 互质。
证毕。
然后推式子:
\begin{aligned}
&(x+1)^n\\
\equiv&(x+1)^{ap+b}\\
\equiv&(x+1)^{ap}\times(x+1)^b\\
\equiv&((x+1)^p)^a\times (x+1)^b\\
\equiv&(x^p+1)^a\times(x+1)^b\\
\equiv&\left(\sum\limits_{i=1}^{a}\binom{a}{i}x^{ip}\right)\times\left(\sum_{j=1}^{b}\binom{b}{j}x^j\right)\pmod p
\end{aligned}
然后考虑如何找到这个式子中的第 m 次项。
设 m=lp+r 。
由于 1\leq j\leq b\leq p ,所以 x^m=x^{lp+r} 只能是由第一个括号里面的 x^{lp} 和第二个括号里面的 x^r 拼凑而来,故 x^m 的系数为 \dbinom{a}{l}\times\dbinom{b}{r} 。
又因为 x^m 本来的系数就是 \dbinom{n}{m} ,所以得:
\dbinom{n}{m}=\dbinom{a}{l}\times\dbinom{b}{r}=\dbinom{\lfloor\frac{n}{p}\rfloor}{\lfloor\frac{m}{p}\rfloor}\times\dbinom{n\bmod p}{m\bmod p}
然后就可以递归求解了。
扩展卢卡斯定理(exLucas)
求:
\binom{n}{m}\bmod p
其中 n,m\leq 10^{18} ,p\leq 10^6 ,不保证 p 为质数 。
先将 p 分解质因数:p=p_1^{k_1}p_2^{k_2}\cdots p_m^{k_m} ,然后求出 \dbinom{n}{m}\bmod p_i^{k_i} ,最后用中国剩余定理合并即可。
求 \dbinom{n}{m}\bmod p^{k} 。
用阶乘拆出来:\dfrac{n!}{m!(n-m)!}\bmod p^k ,注意不能直接求 m! 或 (n-m)! 的逆元,因为不保证它们和 p^k 互质。
设 g(n) 为 n! 中包含的 p 的次数,f(n) 为将 n! 除尽 p 后得到的结果(满足 f(n)=\dfrac{n!}{p^{g(n)}} )。
转化为求:\dfrac{f(n)}{f(m)f(n-m)}\times p^{g(n)-g(m)-g(n-m)}\bmod p^k ,其中 f(n),f(m),f(n-m) 均与 p^k 互质,存在逆元(可以根据扩展欧拉定理 \gcd(a,m)=1\to a^{\varphi(m)}\equiv 1\pmod m 求,或使用 exgcd),于是求出 f 模 p^k 的结果和 g 即可。
求 g(n) 比较简单,不讲了,时间复杂度 O(\log_p n) 。
求 f(n)\bmod p^k 。
对 n! 进行变形:
\begin{aligned}
n!&=1\times 2\times \cdots \times n\\
&=(p\times 2p\times \cdots)(1\times 2\times \cdots )
\end{aligned}
其中左边括号内为所有 p 的倍数,右边括号内为所有与 p 互质的数(当然也与 p^k 互质)。
提出 p 来:
\begin{aligned}
n!&=(p\times 2p\times \cdots)(1\times 2\times \cdots )\\
&=p^{\lfloor\frac{n}{p}\rfloor}\big(\lfloor\frac{n}{p}\rfloor\big)!(1\times 2\times \cdots)
\end{aligned}
那么根据 f(n) 的定义:
f(n)\equiv f\left(\lfloor\frac{n}{p}\rfloor\right)(1\times 2\times \cdots)\pmod{p^k}
现在问题是求出 (1\times 2\times \cdots)\bmod p^k ,为 \leq n 中所有与 p 互质的数的乘积。
求出来之后就可以递归求 f 了,边界条件 f(0)=1 。
求 \left(\prod\limits_{i=1,p\not|i}^ni\right)\bmod p^k 。
法一:
根据 x\equiv x+p^k\pmod{p^k} 可知有长度为 p^k 的循环节,于是等于:
\begin{aligned}
&\left(\prod_{i=1,p\not| i}^{p^k}i\right)^{\lfloor\frac{n}
{p^k}\rfloor}\left(\prod_{i=p^k\lfloor\frac{n}{p^k}\rfloor+1,p\not|i}^{n}i\right)\\
=&\left(\prod_{i=1,p\not| i}^{p^k}i\right)^{\lfloor\frac{n}
{p^k}\rfloor}\left(\prod_{i=1,p\not|i}^{n\bmod p^k}i\right)\\
\end{aligned}
那么可以先 O(p^k) 预处理,每次 O(1) 求解。
法二:
上述方法的瓶颈在求 \leq n 中所有与 p 互质的数的乘积需要 O(p^k) 预处理,这里有另一种适用于 p 小 k 大做法:
设 h_n(x)=\prod\limits_{i=1,p\not|i}^{n}(x+i) ,我们要求的就是 [x^0]h_n(x) 。
可以使用倍增求解:
设 p'=n'p 为 p 的倍数,那么 h_{2p'}(x)=h_{p'}(x)h_{p'}(x+p') 。
考虑归纳证明:假设我们想要知道 [x^i]h_{2p'}(x)\bmod p^{k-i} (0\leq i\leq k ),那么我们只需要知道 [x^i]h_{p'}(x)\bmod p^{k-i} (0\leq i\leq k )。
首先有一个容易证明的引理:对于 A(x)=B(x)C(x) ,假设我们想要知道 [x^i]A(x)\bmod p^{k-i} (0\leq i\leq k ),那么我们只需要知道 [x^i]B(x)\bmod p^{k-i} 和 [x^i]C(x)\bmod p^{k-i} (0\leq i\leq k )。
根据引理,现在只需要证明知道了 [x^i]h_{p'}(x)\bmod p^{k-i} 是否就能知道 [x^i]h_{p'}(x+p')\bmod p^{k-i} 。
设 a_i=[x^i]h_{p'}(x)=q_i\cdot p^{k-i}+r_i 。
\begin{aligned}
\ [x^i]h_{p'}(x+p')&=\sum_{j=i}^k[x^i]a_j(x+p')^j\\
&=\sum_{j=i}^k\binom{j}{i}a_j{p'}^{j-i}\\
&=\sum_{j=i}^k \binom{j}{i}q_jp^{k-i}{n'}^{j-i}+\binom{j}{i}r_j{p'}^{j-i}\\
&\equiv \sum_{j=i}^k \binom{j}{i}r_j{p'}^{j-i}\pmod {p^{k-i}}\\
\end{aligned}
于是要求 [x^i]h_{p'}(x+p')\bmod p^{k-i} 只需 r_j ,即 [x^j]h_{p'}(x)\bmod p^{k-j} ,证毕。
而对于求 h_{2p'+q}(x)=h_{2p'}(x)\prod\limits_{i=2p'+1,p\not|i}^{2p'+q}(x+i) (q<2p )的情况,根据上面提到的引理,我们只需要求出 h_{2p'}(x) 每一位模 p^{k-i} 的多项式,然后暴力乘上后面那部分即可。
具体实现上为了简便,我们可以保留前 k-1 项然后全部对 p^k 取模,而不是对 p^{k-i} 取模。
注意这里求 h_{p'}(x+p') 时需要用到组合数 \dbinom{j}{i} ,这个不能拆成阶乘(因为同样可能不存在逆元),但可以杨辉三角 O(k^2) 预处理。
时间复杂度 O((pk+k^2)\log_pn) 。
如果使用第一种方法:
总时间复杂度 O(\sum p^k)-O(\sum\log_{p^k} n) ,当 p 为质数时取到上界 O(p)-O(\log_p n) 。
特别地,由于我们使用了预处理优化,对于同一个 p 多组询问的复杂度为 O(\sum p^k)-O(T\sum\log_{p^k} n) 。
如果使用第二种方法:
总时间复杂度 O(\sum(pk+k^2)\log_p^2n) 。
Baby Step Giant Step(BSGS,北上广深)
问 a^x\equiv b \pmod p 的最小非负整数解,其中 p\in prime 。
设 t=\sqrt p ,x=it-j 。(因为 p\in prime ,所以 0\leq x\leq p-1 )
\begin{aligned}
a^x&\equiv b \pmod p\\
a^{it-j}&\equiv b\pmod p\\
a^{it}&\equiv b\times a^j \pmod p
\end{aligned}
然后我们先可以 O(\sqrt p) 把 b\times a_j 的所有值存进 hash 表里面。
然后我们再 O(\sqrt p) 枚举 i ,在 hash 表里面找是否有 a^{it}\equiv b\times a^j \pmod p 。
矩阵树定理
对于一个大小为 n^2 的矩阵 A ,定义其行列式 \det(A) 为:(\det(A) 也作 |A| )
\det(A)=\sum\limits_P (-1)^{\mu(P)}\prod\limits_{i=1}^{n}A(i,p_i)
其中 P 是 1\sim n 的一个排列,\mu(P) 是排列 P 的逆序对数。
也就是说,我们在矩阵的每一行选一个数(它们所在的列互不相同),然后乘起来。对于每一种选数方案,赋予它一个与逆序对有关的系数,再加起来。
矩阵树定理有以下的性质:
交换矩阵的两行,行列式变号。
也就等同于交换一个 1\sim n 排列的两个数,整个序列改变的逆序对数为奇数。
证明:
设交换的这两个数位置是在 i 、j 位置,分别是 x 、y 。
我们把序列分成三段,A 段为 1\sim (i-1) ,B 段为 (i+1)\sim (j-1) ,C 段为 (j+1)\sim n 。(如图)
显然,交换 x 和 y 后,x 和 A 、C 两段的相对位置没有改变,y 和 A 、C 两段的相对位置也没有改变,所以由 (A,x) 、(A,y) 、(x,C) 、(y,C) 贡献的逆序对数没有改变。
那么我们只需要考虑 (x,B) 和 (B,y) 两种情况改变的逆序对数就好了。
设 B 段中一共有 S 个数,其中比 x 小的数有 a 个,比 y 大的数有 b 个,那么这段中比 x 大的数有 S-a 个,比 y 小的有 S-b 个。
可以知道交换前,这两种情况贡献的逆序对数共有 a+b 个。
交换后,这两种情况贡献的逆序对数共有 (S-a)+(S-b)=2S-(a+b) 个。
前后贡献的逆序对数之差为 [2S-(a+b)]-(a+b)=2(S-a-b) ,为偶数。
所以 (x,B) 和 (B,y) 两种情况改变的逆序对数为偶数对。
注意到 x 、y 交换后,(i,j) 可能从一开始不是逆序对变成了逆序对,也可能从一开始是逆序对变成了不是逆序对,变动对数为奇数。
所以总的变动对数为奇数对。
证毕。
矩阵的某一行乘上 k ,行列式也乘上 k 。
\begin{vmatrix}
a+a'&b+b'\\
c&d\end{vmatrix}
=\begin{vmatrix}
a&b\\
c&d\end{vmatrix}
+\begin{vmatrix}
a'&b'\\
c&d
\end{vmatrix}
若矩阵中某两行相同,则它的行列式为 0 。
用矩阵的一行加上另一行的倍数,行列式不变。
其中上面的第 2\sim 4 点的证明在这篇题解有详细阐述,不会的可以直接去看。
发现第 5 点操作很像高斯消元中的操作,所以可以考虑用类似高斯消元的方法把行列式变成一个右上三角矩阵,然后这个矩阵的行列式就是对角线上的乘积。
以上是关于行列式和如何求行列式(O(n^3) )。
接下来是矩阵树定理:
无向图。设 A 为邻接矩阵,D 为度数矩阵(D(i,i) 为节点 i 的度数,其他的无值)。则基尔霍夫 (Kirchhoff)矩阵即为 : K=D-A 。把 K 去掉任意第 i 行和第 i 列后的矩阵 K' 的行列式即为答案。
有根扩展。设根为 rt ,则把 K 去掉第 rt 行和第 rt 列后的矩阵 K' 的行列式即为答案。
有向扩展。若为外向树,则度数矩阵统计的是入度;若为内向树,则度数矩阵统计的是出度。
边权扩展。假设有一条边 (u,v,w) ,那么可以看成是在 (u,v) 之间,一共有 w 条重边。
那么我们把度数变成边权就好了。
此时矩阵树定理求的就是 : 所有生成树边权乘积的总和。
积性函数
积性函数:对于一个函数 f(x) ,若有 f(ab)=f(a)f(b) ,其中 \gcd(a,b)=1 ,则称 f(x) 为积性函数。
完全积性函数:对于一个函数 f(x) ,若有 f(ab)=f(a)f(b) ,则称 f(x) 为完全积性函数。
常见的积性函数及其求法:
欧拉函数(\varphi )
定义:\varphi(n) 是小于或等于 n 的正整数中与 n 互质的数的个数。
求法:
首先对于 p\in prime ,肯定有 \varphi(p)=p-1 。
然后对于 p\in prime ,有 \varphi(p^k)=p^{k}-p^{k-1} ,因为除了 p 的倍数以外的其他数都与 p^k 互质。
再加上它是个积性函数,就可以用线性筛求 \varphi 了。
单点求 \varphi
求 \varphi(n) ,设 n=p_1^{k_1}p_2^{k_2}\cdots p_m^{k_m} ,则:
\begin{aligned} \varphi(n)&=\varphi\left(p_1^{k_1}\right) \varphi\left(p_2^{k_2}\right)\cdots\varphi\left(p_m^{k_m}\right)\\ &=p_1^{k_1-1}(p_1-1)\times p_2^{k_2-1}(p_2-1)\times \cdots\times p_m^{k_m-1}(p_m-1)\\ &=\dfrac{n}{p_1p_2\cdots p_m}\times (p_1-1)(p_2-1)\cdots(p_m-1)\\&=n\left(\dfrac{p_1-1}{p_1}\right)\left(\dfrac{p_2-1}{p_2}\right)\cdots\left(\dfrac{p_m-1}{p_m}\right) \end{aligned}
于是就可以 O(\sqrt n) 求 \varphi(n) 了。
莫比乌斯函数(\mu )
\begin{cases}
1&n=1\\
(-1)^k& \text{$n$ 无大于 $1$ 的平方数因数,$n$ 分解质因数为 $n=p_1p_2\cdots p_k$}\\
0 & \text{$n$ 有大于 $1$ 的平方数因数}
\end{cases}
根据这个定义,我们就可以线性筛了。
常见的完全积性函数:
这里记一些可能有用的式子:
证明:对于每个质数幂分别考虑。对于质数 $p$,假设 $i$ 中 $p$ 的幂次为 $a$,$j$ 中 $p$ 的幂次为 $b$。那么左右两式都被统计了 $a+b+1$ 次。
狄利克雷卷积
设 f,g 是两个数论函数,定义他们的狄利克雷卷积是:(f*g)(n)=\sum\limits_{d|n}f(d)g(\dfrac{n}{d}) 。(以下若无特殊说明均用 * 表示狄利克雷卷积)
和普通卷积 c_{i+j}=a_ib_j 不同,这里是 c_{i\times j}=a_ib_j 。
注意狄利克雷卷积针对的是数论函数,而不只是积性函数。
重要的性质:
若 f,g 是积性函数,那么 f*g 也是积性函数。
证明:(其中 \gcd(p,q)=1 )
\begin{aligned}
(f*g)(pq)&=\sum_{d|pq}f(d)g(\frac{pq}{d})\\
&=\sum_{d_1|p}\sum_{d_2|q}f(d_1d_2)g(\frac{pq}{d_1d_2})\\
&=\sum_{d_1|p}f(d_1)g(\frac{p}{d_1})\sum_{d_2|q}f(d_2)g(\frac{q}{d_2})\\
&=(f*g)(p) \times (f*g)(q)
\end{aligned}
其中第一步跳到第二步是因为 p,q 互质,第二步跳到第三步是因为 d_1,d_2 互质且 \dfrac{p}{d_1},\dfrac{q}{d_2} 互质。
若 f,g 是完全积性函数,那么 f*g 也是完全积性函数。
若 f 是积性函数,那么 f^{-1} 也是积性函数。
若 f 是完全积性函数,那么 f^{-1} 也是完全积性函数。
它满足交换律、结合律。
然后对于任意的 f ,有 f*\epsilon=f 。
常见的狄利克雷卷积有:
证明:
若 $n=1$,则易得原式等于 $1$。
若 $n\neq 1$,将 $n$ 质因数分解得 $n=p_1^{a_1}p_2^{a_2}\cdots p_m^{a_m}$。
设枚举到的 $d=p_1^{x_1}p_2^{x_2}\cdots p_m^{x_m}
按枚举到的 d 分类讨论:
对于所有的 x_i ,存在有至少一个 x_i>1 。
此时 \mu(d)=0 。
对于所有的 x_i ,都有 x_i\leq 1 。
那么这部分的 ans 为:
ans=\sum_{i=0}^{m}\binom{m}{i}(-1)^i
其中第一个 \sum\limits_{i=0}^m 枚举所有的 x_j 中(1\leq j \leq m ),有 i 个 x_j 大于 0 ,即有 i 个 x_j 为 1 。
然后 \dbinom{m}{i} 意思是从所有的 x_j 中(1\leq j \leq m ),选出 i 个 x_j 为 1 的方案数。
然后 (-1)^i 就是 \mu(d) 。
至于 ans 为什么等于 0 :
若 m\equiv 1\pmod 2 。
则 +\dbinom{m}{0} 会和 -\dbinom{m}{m} 抵消,-\dbinom{m}{1} 会和 +\dbinom{m}{m-1} 抵消……以此类推,可得 ans=0 。
若 m\equiv 0 \pmod 2 。
将每个 \dbinom{m}{i} 都画在杨辉三角里:
\begin{aligned}ans=&\binom{m}{0}-\binom{m}{1}+\binom{m}{2}-\cdots-\binom{m}{m-1}+\binom{m}{m}\\=&\binom{m-1}{0}-\left[\binom{m-1}{0}+\binom{m-1}{1}\right]+\left[\binom{m-1}{1}+\binom{m-1}{2}\right]-\cdots\\&-\left[\binom{m-1}{m-2}+\binom{m-1}{m-1}\right]+\dbinom{m-1}{m-1}\\=&0\end{aligned}
所以无论如何,都有 ans=0 。
所以当 n\neq 1 时,原式为 0 。
证毕。
证明:
首先把式子转换一下:$\sum\limits_{d|n}\varphi(d)=\sum\limits_{d|n}\varphi(\dfrac{n}{d})$。
考虑每一个 $\varphi(\dfrac{n}{d})$ 代表着什么。
可以把 $\varphi(\dfrac{n}{d})$ 看成小于等于 $n$ 的所有正整数中,与 $n$ 的 $\gcd$ 为 $d$ 的个数。
那么 $\sum\limits_{d|n}\varphi(\dfrac{n}{d})$ 就可以看成所有小于等于 $n$ 的正整数个数,得 $\sum\limits_{d|n}\varphi(\dfrac{n}{d})=n$。
证毕。
证明:
首先由(2)得 $\varphi*I=id$。
则 $\varphi*I*\mu=id*\mu$。
则 $\varphi*\epsilon=id*\mu$。
即 $\varphi=id*\mu$。
证毕。
贝尔级数
定义数论函数 f 在质数 p 上的贝尔级数 f_p(x) :
f_p(x)=\sum_{i\geq 0}f(p^i)x^i
那么有:
(f*g)_p(x)=f_p(x)g_p(x)
证明:
\begin{aligned}
(f*g)_p(x)=&\sum_{i\geq 0}(f*g)(p^i)x^i\\
=&\sum_{i\geq 0}x^i\sum_{d|p^i}f(d)g(\frac{p^i}{d})\\
=&\sum_{i\geq 0}\sum_{j=0}^if(p^{j})x^jg(p^{i-j})x^{i-j}\\
=&f_p(x)g_p(x)
\end{aligned}
同时,若 f_p(x)g_p(x)=h_p(x) 且 h 为积性函数,那么 h=f*g 。
这样就能用普通的卷积来刻画狄利克雷卷积。
因为狄利克雷卷积 (f*g)(n) 的实际上是在枚举两个数 a,b 使得 ab=n ,那么相当于对于 n 的每一个质因子及其幂 p^k 枚举 i+j=k ,那就由狄利克雷卷积变成了普通卷积。
对于完全积性函数 f 有:
f_p(x)=\sum_{i\geq 0}f(p^i)x^i=\sum_{i\geq 0}f^i(p)x^i=\dfrac{1}{1-f(p)x}
常见数论函数的贝尔级数:
于是你可以看出这些数论函数之间的很多关系,如:
通过贝尔级数也可以更加容易地看出一些奇怪的数论函数与常见的数论函数的关系,这对杜教筛的构造有很大的作用。
莫比乌斯反演
定理:若 g=f*I ,即 g(n)=\sum\limits_{d|n}f(d) ,则 f=g*\mu ,即 f(n)=\sum\limits_{d|n}g(d)\mu(\dfrac{n}{d}) 。
证明: g=f*I 得 g*\mu=f*I*\mu=f*\epsilon=f ,证毕。
杜教筛
求一个积性函数 f 的前缀和,记 S(n)=\sum\limits_{i=1}^n f(i) 。
再找一个积性函数 g ,先考虑它们的狄利克雷卷积的前缀和:
\begin{aligned}
&\sum_{i=1}^{n}(f*g)(i)\\
=&\sum_{i=1}^n\sum_{d|i}f(d)g(\frac{i}{d})\\
=&\sum_{d=1}^ng(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i)\\
=&\sum_{d=1}^ng(d)S\left(\lfloor\frac{n}{d}\rfloor\right)
\end{aligned}
考虑 g(1)S(n) 是什么,发现:
\begin{aligned}
g(1)S(n)&=\sum_{d=1}^ng(d)S\left(\lfloor\frac{n}{d}\rfloor\right)-\sum_{d=2}^ng(d)S\left(\lfloor\frac{n}{d}\rfloor\right)\\
&=\sum_{i=1}^{n}(f*g)(i)-\sum_{d=2}^ng(d)S\left(\lfloor\frac{n}{d}\rfloor\right)
\end{aligned}
那么我们就可以找到或构造一个 g ,使得 \sum\limits_{i=1}^{n}(f*g)(i) 和 g 的前缀和可以快速算,然后数论分块递归求解。
然后可以预处理出前若干个答案的前缀和,并且加上记忆化,可以把时间复杂度优化成 O(n^{\frac{2}{3}}) ,其中要处理前 n^{\frac{2}{3}} 个答案的前缀和。
min_25 筛
记号:
记质数集合为 \mathbb{P} 。
如无特殊说明,以下记为 p 的变量的取值集合为质数集合。
为了方便,有时用 (a/b) 表示 \left\lfloor\dfrac{a}{b}\right\rfloor 。
特别地,p_0=1 ,\operatorname{lpf}(1)=1 ,那么就保证了不存在 i 使得 \operatorname{lpf}(1)>p_i 。
min25 筛用来解决这样一类问题:问 f(x) 前缀和 $\sum\limits {i=1}^n f(i)。其中 f(x)$ 需要满足以下条件:
具体为什么要满足这些条件在下面有体现。
我们将 i=1\sim n 分为质数和合数两类讨论:
\begin{aligned}
\sum_{i=1}^nf(i)=f(1)+\sum_{p\leq n}f(p)+\sum_{p\leq \sqrt n,p^e\leq n}f(p^e)\left(\sum_{2\leq i\leq (n/p^e),\operatorname{lpf}(i)>p}f(i)+[e>1]\right)
\end{aligned}
注意若 x 为合数,那么 x 的最小质因数 p\leq \sqrt x ,所以我们 p 只需枚举到 \sqrt n 即可。
设 S(n,i)=\sum\limits_{x=1}^n[\operatorname{lpf}(x)>p_i]f(x) ,那么就可以递推:
S(n,i)=\sum_{p\leq n}f(p)-\sum_{p\leq p_i}f(p)+\sum_{k>i,p_k\leq \sqrt n,p_k^e\leq n}f(p_k^e)\left(S\big(\lfloor\frac{n}{p_k^e}\rfloor,k\big)+[e>1]\right)
根据定义式,f(1)+S(n,0) 就是答案。
如果记 F_{\operatorname{prime}}(n)=\sum\limits_{p\leq n}f(p) ,那么:
S(n,i)=F_{\operatorname{prime}}(n)-F_{\operatorname{prime}}(p_i)+\sum_{k>i,p_k\leq \sqrt n,p_k^e\leq n}f(p_k^e)\left(S\big(\lfloor\frac{n}{p_k^e}\rfloor,k\big)+[e>1]\right)
我们需要解决 F_{\operatorname{prime}}(n) 的问题:
假设 f(p) 为一多项式,设 f(p)=\sum\limits_{k=0}^ma_kp^k 。
F_{\operatorname{prime}}(n)=\sum_{p\leq n}f(p)=\sum_{p\leq n}\sum_{k=0}^ma_kp^k=\sum_{k=0}^ma_k\sum_{p\leq n}p^k
于是我们只需要求出 \sum\limits_{p\leq n}p^k 即可。
定义函数 g_k(n,i) 为:
g_k(n,i)=\sum_{x=1}^n[(x\in \mathbb{P})\lor (\operatorname{lpf}(x)>p_i)]x^k
尝试递推 g_k(n,i) :
\begin{aligned}
g_k(n,i)&=\sum_{x=1}^n[(x\in \mathbb{P})\lor (\operatorname{lpf}(x)>p_i)]x^k\\
&=\sum_{x=1}^n[(x\in \mathbb{P})\lor (\operatorname{lpf}(x)>p_{i-1})]x^k-\sum_{x=1}^n[(x\not \in \mathbb{P})\land (\operatorname{lpf}(x)=p_i)]x^k\\
&=g_k(n,i-1)-\sum_{y=1}^{(n/p_i)}[\operatorname{lpf}(y)\ge p_i](yp_i)^k\\
&=g_k(n,i-1)-p_i^k\sum_{y=1}^{(n/p_i)}[\operatorname{lpf}(y)> p_{i-1}]y^k\\
&=g_k(n,i-1)-p_i^k\left(\sum_{y=1}^{(n/p_i)}[(y\in \mathbb{P})\lor(\operatorname{lpf}(y)> p_{i-1})]y^k-\sum_{y=1}^{(n/p_i)}[(y\in \mathbb{P})\land(\operatorname{lpf}(y)\leq p_{i-1})]y^k\right)\\
&=g_k(n,i-1)-p_i^k\left(\sum_{y=1}^{(n/p_i)}[(y\in \mathbb{P})\lor(\operatorname{lpf}(y)> p_{i-1})]y^k-\sum_{y=1}^{p_{i-1}}[y\in \mathbb{P}]y^k\right)\\
&=g_k(n,i-1)-p_i^k\left(\sum_{y=1}^{(n/p_i)}[(y\in \mathbb{P})\lor(\operatorname{lpf}(y)> p_{i-1})]y^k-\sum_{y=1}^{p_{i-1}}[(y\in \mathbb{P})\lor(\operatorname{lpf}(y)> p_{i-1})]y^k\right)\\
&=g_k(n,i-1)-p_i^k\bigg(g_k\big((n/p_i),i-1\big)-g_k\big(p_{i-1},i-1\big)\bigg)\\
\end{aligned}
大概意思就是说,g_k(n,i) 一定是由 g_k(n,i-1) 减去最小质因数恰好是 p_i 的合数 x 的贡献得到,于是枚举最小质因数大于等于 p_i 的数 y 使得 x=p_iy ,也就是 g_k((n/p_i),i-1) ,但还需减去 y 为小于 p_i 的质数的贡献,即 g_k(p_{i-1},i-1) 。
注意倒数第 4 行推到倒数第 3 行时是默认了 p_{i-1}\leq (n/p_i)=\lfloor\dfrac{n}{p_i}\rfloor 的,因为如果 p_{i-1}>\lfloor\dfrac{n}{p_i}\rfloor 的话,就有 p_{i-1}+1>\lfloor\dfrac{n}{p_i}\rfloor+1>\dfrac{n}{p_i} ,于是 \dfrac{n}{p_i}<p_{i-1}+1<p_i ,则 p_i^2>n ,就不存在最小质因数恰好是 p_i 的合数 x ,此时直接 g_k(n,i)=g_k(n,i-1) 即可。
那么:
g_k(n,i)=
\begin{cases}
g_k(n,i-1)&\text{if }p_i^2>n\\
g_k(n,i-1)-p_i^k\bigg(g_k\big((n/p_i),i-1\big)-g_k\big(p_{i-1},i-1\big)\bigg)&\text{if }p_i^2\leq n
\end{cases}
于是只要我们知道 g_k(n,0) 就可以递推了。
g_k(n,0)=\sum_{x=1}^n[(x\in \mathbb{P})\lor (\operatorname{lpf}(x)>p_0)]x^k=\sum_{x=2}^nx^k
这是自然数幂求和,直接 O(k) 做即可。
于是 \sum\limits_{p\leq n}p^k=g_k(n,\pi(\sqrt n)) 即可。(其中 \pi(n) 表示小于等于 n 的质数个数)
于是 F_{\operatorname{prime}}(n) 就被我们解决了,于是就可以递推 S(n,i) 了。
时间复杂度不会……
口胡 min_25 筛的主要思路(主要方便自己记):
求 \sum\limits_{i=1}^nf(i) ,我们先对 1\sim n 分成质数和合数,合数的部分我们通过枚举最小质因数递推,关键在于质数的部分,即 \sum\limits_{p\leq n}f(p) 。然后 f(p) 是低次多项式,可以拆成不同次分开处理,即求 \sum\limits_{p\leq n}p^k 。然后发现这个时候 \sum\limits_{i=1}^ni^k 是好算的,于是又反过来,用全体数减去合数得到质数的答案,合数的部分也是枚举最小质因数递推。
一些代码实现的细节:
注意到求 S(x,y) 时,向下递归的形式是 S((x/b),\cdots) ,其中 b 为某一个正整数。于是我们就猜想:是否所有求 S(x,y) 时,x 都能表示为 \left\lfloor\dfrac{n}{c}\right\rfloor 的形式。答案当然是能。g(x,y) 也是同理。
于是出现在 S 和 g 第一维的 x 的取值只能是:(n/1),(n/2),(n/3),\cdots,(n/n) ,共 O(\sqrt n) 种。
那么我们整除分块,将这些取值离散化,然后就可以把空间复杂度降到只有 O(\sqrt n) 。将每一种取值给一个标号,记取值 v 的标号为 id(v) 。
显然如果我们真的要存 id(v) 的话就需要用 map,但我们有更好的方法。
具体来说,对于 x=(n/i) ,若 x\leq \sqrt n ,那么我们令 id_1(x)\gets id(x) ;若 x>\sqrt n ,那么我们令 id_2(n/x)\gets id(x) 。
查询 x 的标号的时候我们就直接判断 x 与 \sqrt n 的大小关系再进入 id_1 或 id_2 查询即可。
多项式
下面一部分专题都会与多项式相关。
原根
参考了 cmd 的博客。
不得不说 cmd 是真的卷王。
阶 :对于任意的有限群 G=\{e,a_1,\cdots,a_n\} ,对于任意的 a\in G 都存在一个正整数常数 \operatorname{ord}(a) 使得 a^{\operatorname{ord}(a)}=e ,称 \operatorname{ord}(a) 为 a 的阶。
证明:
引理:对于群 G 以及 a,b\in G ,若 a\cdot b=a ,那么 b 是 G 的单位元。
证明:只需证对于其他的 a'\in G 也满足 a'\cdot b=a' 即可。
构造 a,a^2,\cdots,a^{n+1} ,由于这 n+1 个元素都属于 G ,所以必有二者相等,设为 a^x=a^y(x<y) ,那么 a^x\cdot a^{y-x}=a^x 。再根据引理可证 a^{y-x} 为单位元。
意思是 a 的若干次幂在有限群 G 中一定至少构成一个 \rho 字形,那么 a 的环的大小次幂就是单位元。
推论 1:a 的若干次幂不仅是 \rho 字形,而且就是一个环,\operatorname{ord}(a) 就是环的大小。
推论 2:a,a^2,\cdots,a^{\operatorname{ord}(a)} 两两不同。
证明:由上面证明过程容易得到。
整数模 n 乘法群(模 n 既约剩余类) :
这里我们用 [a] 表示模 n 与 a 同余的剩余类。
如果一个模 n 的剩余类 [a] 中的所有数都与 n 互质(等价于 (n,a)=1 ),那么我们称 [a] 为一个模 n 既约剩余类。
所有的模 n 既约剩余类在剩余类的乘法运算下组成一个乘法群,称为整数模 n 乘法群,姑且记作 G_n 。
容易验证 G_n 满足群公理:
封闭性:若 [a],[b]\in G_n ,则 ab 也与 n 互质,故 [a][b]=[ab]\in G_n 。
结合律:由乘法运算法则不难得到。
存在单位元:[1] 肯定属于 G_n ,它是单位元。
存在逆元:即对于任意 [a]\in G_n ,都存在 (x,n)=1 且 ax\equiv 1\pmod n 。
转化为 ax+ny=1 ,由裴蜀定理可知有解当且仅当 (a,n)=1 ,得证。
下面为了方便和直观,建议你可以直接将模 n 余 a 剩余类当成数字 a ,即直接将 G_n 看成所有小于等于 n 且与 n 互质的正整数在模 n 意义的乘法下组成的群。
那么对于 a\in G_n ,它的阶就是最小的正整数 \operatorname{ord}(a) 使得 a^{\operatorname{ord}(a)}\equiv 1\pmod p ,即幂的最小循环节,由群的封闭性可知它的上界为 |G_n|=\varphi (n) 。
循环群 :若一个群 G 的每一个元都是 G 的某一个固定元 a 的若干次幂,则称 G 为循环群,记作 G=\{a^m|m\in \mathbb{Z}\} 。a 称为 G 的一个生成元。
容易看出 \operatorname{ord}(a)=|G| ,也容易得到一个群是循环群当且仅当存在 a\in G 满足 \operatorname{ord}(a)=|G| 。
循环群分两种:
若 |G|=\infty ,则 G 的生成元有且仅有两个,分别为 a 和 a^{-1} 。
证明:此时 G 与整数加法群同构(考虑幂乘法和指数加法的对应关系),而整数加法群只有两个生成元 1 和 -1 。
若 |G|=n ,则 G 的生成元有 \varphi(n) 个,它们为 \{a^k|(k,n)=1\} 。
证明:考虑 a 的幂,这构成一个大小为 n 的、恰好包含 G 中所有元素的环。
那么 a^k 的若干次幂相当于在这个环上走若干次每次走 k 步,容易知道这样走也会构成一个大小为 \frac{n}{\gcd(k,n)} 的环。
所以若 (k,n)=1 ,那么 a^k 一定是生成元。否则 a^k 一定不是生成元。
而且 G 中的任意一个元素都可以被表示成 a^k 的形式,所以所有可能的生成元都被我们考虑到了。
推论:设 b\in G 且 b=a^k ,那么 \operatorname{ord}(b)=\frac{n}{\gcd(k,n)} 。证明由上面证明过程易得。
接下来需要考虑的是第二种循环群。
那么若 p 有原根,个数为 \varphi(|G_p|)=\varphi(\varphi(p)) 。若一个数与 p 互质,它肯定能用原根的若干次幂表示。
一个结论是当 p=2,4,p^c,2p^c (p 为奇质数)时,G_p 才是循环群,p 才有原根。
接下来的问题是:假设我们已经判定了 p 有原根,然后如何求 p 的原根。
一个结论是 p 的最小原根不超过 p^{0.25} 。
那么我们可以从小往大枚举数 g 满足 (g,p)=1 ,并判断它是否为 p 的原根,只需判定它的阶是否为 \varphi(p) 即可。
枚举阶快速幂判定肯定是慢的,更优美的做法是:我们并不需要求出 g 的阶的准确值,只需要知道它的阶是否为 \varphi(p) 即可。
由循环群中的推论可知,g 的阶一定是 \varphi(p) 的因数。
又注意到若 g^k\equiv 1\pmod p ,那么 g^{ik}\equiv 1\pmod p 。
所以 g 的阶为 \varphi(p) 的严格因数的判定条件就可以转化为 \exist i\in[1,m],g^{\tfrac{p}{p_i}}\equiv 1\pmod p ,其中 p_1,\cdots,p_m 为 p 的质因数。
于是单次判定某个 g 的复杂度就降为 O(\omega(p)\log p) ,找某个原根的总复杂度就是 O(p^{0.25}\omega(p)\log p) 。
Bluestein 算法
给出多项式 A(x)=\sum\limits_{i=0}^na_ix^i 和 c ,求 A(1),A(c),A(c^2),\cdots 的算法。
FFT 做的是 c 为 2 的幂次的单位根的特殊情况。而这种算法对 c 没有特殊性质。
关键在于 ki=\dbinom{k+i}{2}-\dbinom{k}{2}-\dbinom{i}{2} ,于是:
\begin{aligned}
A(c^k)&=\sum_{i=0}^n a_ic^{ki}\\
&=\sum_{i=0}^na_ic^{\tbinom{k+i}{2}-\tbinom{k}{2}-\tbinom{i}{2}}\\
&=c^{-\tbinom{k}{2}}\sum_{i=0}^na_ic^{-\tbinom{i}{2}}\cdot c^{\tbinom{k+i}{2}}
\end{aligned}
减法卷积即可。
多项式求逆
给定多项式 F(x) ,求多项式 G(x) 使得 F(x)G(x)\equiv1\pmod {x^n} ,其中 F(x)G(x) 的系数对一个数 p 取模。(对于一个多项式 F(x) ,F(x)\bmod x^n 的意思是只看 F(x) 中的第 0 项到第 n-1 项)
用递归求解。(以下的所有多项式都是在系数模了 p 意义下的)
首先假设我们已经求出了 F(x)H(x)\equiv 1 \pmod{x^{\left\lceil\frac{n}{2}\right\rceil}} (记为 1 式),现在要求 F(x)G(x)\equiv 1 \pmod{x^n} 。
首先如果有 F(x)G(x)\equiv 1 \pmod{x^n} ,那么一定有 F(x)G(x)\equiv 1 \pmod{x^{\left\lceil\frac{n}{2}\right\rceil}} (记为 2 式),因为等号右边的那个 1 只能在常数项出现,其他项的系数都是 0 。
由 1 式减 2 式得 F(x)\left(H(x)-G(x)\right)\equiv 0\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}} 。
由于 F(x)\not\equiv 0\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}} ,所以 H(x)-G(x)\equiv 0 \pmod{x^{\left\lceil\frac{n}{2}\right\rceil}} 。
然后等式两边平方,因为 H(x)-G(x)\equiv 0 \pmod{x^{\frac{n}{2}}} ,即 H(x)-G(x) 的第 0 项到第 \left\lceil\dfrac{n}{2}\right\rceil-1 项系数都是 0 ,那么 (H(x)-G(x))^2 肯定满足第 0 项到第 n-1 项系数都是 0 ,所以有:
H(x)^2+G(x)^2-2G(x)H(x)\equiv 0\pmod{x^n}
两边同时乘 F(x) 得:
\begin{aligned}
F(x)H(x)^2+F(x)G(x)^2-2F(x)G(x)H(x)&\equiv 0\pmod{x^n}\\
F(x)H(x)^2+G(x)-2H(x)&\equiv 0\pmod{x^n}\\
\end{aligned}
上面第二步的化简是用 F(x)G(x)\equiv 1 \pmod{x^n} 化简的,注意不能用 F(x)H(x)\equiv 1 \pmod{x^{\left\lceil\frac{n}{2}\right\rceil}} ,因为两个式子 \bmod 的东西不一样。
那么就可以得到:G(x)=2H(x)-F(x)H(x)^2 。
那么我们就能通过 \left\lceil\dfrac{n}{2}\right\rceil 的答案求出 n 时的答案了,所以一直向下递归。
最后当递归到 1 的时候,要求 F(x)G(x)\equiv 1\pmod{x^1} ,也就是 F(x) 的常数项和 G(x) 的常数项乘起来模 p 为 1 ,那么只用求 F(x) 常数项的逆元就好了。
总时间复杂度为:(设 n=2^m ,以下的 \log 均以 2 为底)
\sum\limits_{i=0}^{m-1}2^i\log(2^i)<\sum\limits_{i=0}^{m-1}2^i\log(n)=\left(2^{m-1+1}-1\right)\log(n)=(n-1)\log(n)<n\log n
多项式开根
跟多项式求逆差不多,用递归求解。
假设我们现在要求 G(x) 使得 G(x)^2\equiv F(x)\pmod{x^n} ,已经求出了 H(x)^2\equiv F(x)\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}} 。(记为 1 式)
显然若 G(x)^2\equiv F(x)\pmod{x^n} ,则 G(x)\equiv \sqrt{F(x)}\equiv H(x)\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}} 。
\begin{aligned}
G(x)&\equiv H(x)\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\\
G(x)-H(x)&\equiv 0\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\\
(G(x)-H(x))^2&\equiv 0\pmod{x^n}\\
G(x)^2+H(x)^2-2G(x)H(x)&\equiv 0\pmod{x^n}\\
F(x)+H(x)^2-2G(x)H(x)&\equiv 0\pmod{x^n}\\
G(x)&\equiv\frac{F(x)+H(x)^2}{2H(x)}\pmod{x^n}
\end{aligned}
然后就能递归求逆做了。
时间复杂度 O(n\log n) 。
多项式求导
设 f(x)=\sum\limits_{i\geq0}a_ix^i ,则:
\begin{aligned}
f'(x)=&\lim_{d\rightarrow 0}\sum_{i\geq0}\frac{a_i(x+d)^i-a_ix^i}{d}\\
=&\lim_{d\rightarrow 0}\sum_{i\geq0}a_i\frac{(x+d)^i-x^i}{d}
\end{aligned}
我们现在只看 \dfrac{(x+d)^i-x^i}{d} 怎么求:
\begin{aligned}
&\frac{(x+d)^i-x^i}{d}\\
=&\frac{\sum\limits_{j=0}^i\dbinom{i}{j}x^jd^{i-j}-x^i}{d}\\
=&\frac{\sum\limits_{j=0}^{i-1}\dbinom{i}{j}x^jd^{i-j}+x^id^0-x^i}{d}\\
=&\frac{\sum\limits_{j=0}^{i-1}\dbinom{i}{j}x^jd^{i-j}}{d}\\
=&\sum\limits_{j=0}^{i-1}\dbinom{i}{j}x^jd^{i-j-1}
\end{aligned}
则 j<i-1 时,i-j-1>0 ,又由于 d 趋近于 0 ,则 d^{i-j-1}=0 。
所以:
\sum\limits_{j=0}^{i-1}\dbinom{i}{j}x^jd^{i-j-1}=\dbinom{i}{i-1}x^{i-1}d^0=ix^{i-1}
那么:
\begin{aligned}
f'(x)=&\lim_{d\rightarrow 0}\sum_{i\geq0}a_i\frac{(x+d)^i-x^i}{d}\\
=&\sum_{i\geq 1}a_iix^{i-1}
\end{aligned}
同样的,若知道 f'(x)=\sum\limits_{i\geq 0}a_ix^i ,那么 f(x)=\sum\limits_{i\geq 1}\dfrac{a_{i-1}}{i}x^i 。
多项式 \ln
设 F(x)=\ln x ,那么我们要求的即为 G(x)\equiv F(A(x))\pmod{x^n} 。
对两边同时求导:
G'(x)\equiv F'(A(x))A'(x)\equiv\dfrac{A'(x)}{A(x)}\pmod{x^n}
直接求出 G'(x) 再积分即可。
多项式 \exp
泰勒展开
我们为了近似表达一函数 f(x) ,构造另一函数 g(x) ,使得 f(x)=g(x),f'(x)=g'(x),f^{(n)}(x)=g^{(n)}(x) 。
设 g(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n 。
由于 f(x)=g(x) ,代入 x=0 可得 f(0)=g(0) ,又 g(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n ,得 a_0=f(0) ;
由于 f'(x)=g'(x) ,代入 x=0 可得 f'(0)=g'(0) ,又 g'(x)=1a_1+2a_2x+\cdots+na_nx^{n-1} ,得 a_1=\dfrac{f'(0)}{1!} ;
由于 f''(x)=g''(x) ,代入 x=0 可得 f''(0)=g''(0) ,又 g''(x)=2a_2+6a_3x\cdots+n(n-1)a_nx^{n-2} ,得 a_2=\dfrac{f''(0)}{2!} ;
……
最后得到 a_n=\dfrac{f^{(n)}(0)}{n!} ,故:
g(x)=f(0)+\dfrac{f'(0)}{1!}x+\dfrac{f''(0)}{2!}x^2+\cdots+\dfrac{f^{(n)}(0)}{n!}x^n
这是代入 x=0 的特殊情况,而我们代入 x=x_0 也是同理的,有:
g(x)=f(x_0)+\dfrac{f'(x_0)}{1!}(x-x_0)+\dfrac{f''(x_0)}{2!}(x-x_0)^2+\cdots+\dfrac{f^{(n)}(x_0)}{n!}(x-x_0)^n
当 n\to \infty 时有 f(x)=g(x) 。
多项式牛顿迭代
求一类方程的解:
F(G(x))\equiv 0 \pmod{x^n}
假设我们已经得知了 G_0(x)\equiv G(x)\pmod {x^{\lceil\frac{n}{2}\rceil}} ,要求 G(x)\bmod x^n 。
将 F(G(x)) 在 G_0(x) 处泰勒展开:
F(G(x))=\sum_{n=0}^{\infty}\dfrac{F^{(n)}(G_0(x))}{n!}(G(x)-G_0(x))^n
代入 G(x)=G_1(x) 得:
F(G_1(x))=\sum_{n=0}^{\infty}\dfrac{F^{(n)}(G_0(x))}{n!}(G_1(x)-G_0(x))^n\equiv 0\pmod{x^n}
又由于 G_0(x)\equiv G(x)\pmod {x^{\lceil\frac{n}{2}\rceil}} ,则:
\begin{aligned}
G_0(x)&\equiv G(x)\pmod {x^{\lceil\frac{n}{2}\rceil}}\\
G_0(x)-G(x)&\equiv 0\pmod {x^{\lceil\frac{n}{2}\rceil}}\\
(G_0(x)-G(x))^2&\equiv 0\pmod {x^{n}}\\
\end{aligned}
所以次数大于等于 2 的项都被消掉了。
那么:
\begin{aligned}
\sum_{n=0}^{\infty}\dfrac{F^{(n)}(G_0(x))}{n!}(G(x)-G_0(x))^n&\equiv 0\pmod{x^n}\\
F(G_0(x))+F'(G_0(x))(G(x)-G_0(x))&\equiv 0\pmod{x^n}\\
\end{aligned}
当 F'(G_0(x)) 常数项都非零时,得到迭代公式:
G(x)\equiv G_0(x)-\dfrac{F(G_0(x))}{F'(G_0(x))}\pmod{x^n}
当 F'(G_0(x)) 末 k 项为 0 且第 k+1 项非零时(此时 F(G_0(x)) 末 k 项也应为 0 ),一般来说 k 都是常数,那么牛顿迭代的时候考虑由 G(x)\bmod x^{k+\lceil\frac{n}{2}\rceil} 推到 G(x)\bmod x^{k+n} 。
这是因为对于 A(x)B(x)\equiv C(x)\pmod{x^n} ,其中 A,C 末 k 项为 0 且第 k+1 项非零时,只能得到 B(x)\equiv\frac{C(x)}{A(x)}\pmod{x^{n-k}} 。
多项式 \exp
\begin{aligned}
B(x)&\equiv e^{A(x)}\pmod {x^n}\\
\ln B(x)-A(x)&\equiv 0\pmod {x^n}
\end{aligned}
那么就是要求函数 F(G(x))=\ln G(x)-A(x) 的零点,即 F(G(x))\equiv 0\pmod{x^n} 的 G(x) 的解。
套用牛顿迭代,假设我们已经得知了 G_0(x)\equiv G(x)\pmod {x^{\lceil\frac{n}{2}\rceil}} ,要求 G_1(x)\equiv G(x)\pmod {x^n} 。
根据迭代公式,得:
\begin{aligned}
G_1(x)&\equiv G_0(x)-\dfrac{F(G_0(x))}{F'(G_0(x))}\pmod{x^n}\\
&\equiv G_0(x)-\dfrac{\ln G_0(x)-A(x)}{\dfrac{1}{G_0(x)}}\pmod{x^n}\\
&\equiv G_0(x)(1-\ln G_0(x)+A(x))\pmod{x^n}
\end{aligned}
那么就可以求了。
时间复杂度 O(n\log n) 。
多项式带余除法
给定一个 n 次多项式 F(x) 和一个 m 次多项式 G(x) ,要求求出多项式 Q(x) 和 R(x) ,满足下面的条件:
考虑先求出商,再求余数。
设 A(x) 为一 n 次多项式,定义 A_r(x)=x^nA(\dfrac{1}{x}) ,易知 A_r[i]=A[n-i] 。
\begin{aligned}
F(x) &= Q(x) \times G(x)+ R(x)\\
F(\frac{1}{x}) &= Q(\frac{1}{x}) \times G(\frac{1}{x})+ R(\frac{1}{x})\\
x^nF(\frac{1}{x}) &= x^{n-m}Q(\frac{1}{x}) \times x^{m}G(\frac{1}{x})+ x^nR(\frac{1}{x})\\
F_r(x)&=Q_r(x)\times G_r(x)+x^{n-m+1}R_r(x)\\
F_r(x)&\equiv Q_r(x)\times G_r(x)\pmod{x^{n-m+1}}\\
Q_r(x)&\equiv G_r(x)\times F_r^{-1}(x)\pmod{x^{n-m+1}}
\end{aligned}
然后直接求出 Q_r(x) ,再求出 Q(x) 和 R(x) 即可。
时间复杂度 O(n\log n) 。
多项式多点求值
给定一个 n 次多项式 F(x) ,现在请你对于 i \in [1,m] ,求出 F(a_i) 。
注意到多项式取模 F(x)=Q(x)G(x)+R(x) 中,若 G(x_0)=0 ,则 R(x_0)=F(x_0) 。
所以我们考虑对于 i 构造 G_i(x)=x-a_i ,那么显然 F(a_i)=F(a_i)\bmod G_i(a_i) 。
所以我们只需要对于每一个 i 求出 F(x)\bmod G_i(x) 即可,注意到此时模出来是一个常数,那么就是 F(a_i) 。不妨设 F(x)=Q_i(x)G_i(x)+R_i(x) 。
考虑如何求 R_i(x) 。
考虑分治,设当前要求出对于所有 i\in [L,R] 的 F(x)\bmod G_i(x) ,那么我们先让 F(x) 对 \prod\limits_{i=L}^RG_i(x) 取模,然后在往 [L,mid] 、[mid+1,R] 递归。容易得知这样做不会影响答案,而且每次能将 F(x) 的次数减少一半,时间复杂度 O(n\log n+m\log^2 m) (带上 n\log n 的原因是一开始还要将 F(x) 对 \prod\limits_{i=1}^mG_i(x) 取模,这一次消耗的时间是 O(n\log n) 的)。
但鉴于多项式取模的极大常数,这个可能连 10^5 都跑不过去。
有一种常数更小的做法:
重新分析一下,R_i(x)=F(x)-Q_i(x)G_i(x) ,而 R_i(x) 只有常数项,所以如果我们知道了 Q_i(x) 的常数项,就能知道 R_i(x) 了。
考虑多项式取模是如何求出 Q_i(x) 的:Q_{i,r}(x)\equiv F_r(x)G_{i,r}^{-1}(x)\pmod{x^{n-m+1}} ,即 Q_{i,r}(x)\equiv F_r(x)G_{i,r}^{-1}(x)\pmod{x^{n}} (这里下标带 r 是指反转)。
仍然考虑分治,设 G_s(x)=\prod\limits_{i=L}^R G_i(x) ,G_0(x)=\prod\limits_{i=L}^{mid} G_i(x) ,G_1(x)=\prod\limits_{i=mid+1}^R G_i(x) 。
那么易知 G_s(x)=G_0(x)G_1(x) ,则 G_s(\dfrac{1}{x})=G_0(\dfrac{1}{x})G_1(\dfrac{1}{x}) 、x^{R-L+1}G_s(\dfrac{1}{x})=x^{mid-L+1}G_0(\dfrac{1}{x})x^{R-mid}G_1(\dfrac{1}{x}) ,即 G_{s,r}(x)=G_{0,r}(x)G_{1,r}(x) 。
那么假设我们已经知道了 F_r(x)G_{s,r}^{-1}(x) 的结果,那么我们只需要对这个结果乘上 G_{1,r}(x) 即可知道 F_r(x)G_{0,r}^{-1}(x) 的值;对这个结果乘上 G_{0,r}(x) 即可知道 F_r(x)G_{1,r}^{-1}(x) 的值。递归到底层即可。
但这样时间复杂度是不对的,因为我们为了保证底层结果的正确,每次乘法都要在模 x^n 意义下进行。
但注意到我们只需要求 Q_{i}(x) 的常数项,即 [x^{n-1}]Q_{i,r}(x) ,也就是说只有某一项可能对最高项产生贡献才需要保留,否则可以直接舍弃。
那么假设从当前到底层还要乘上的 G(x) 之积的最高次为 k 次,那么当前的 Q(x) 只需要保留最高的 k+1 项,因为剩下的项都不会对最后的 Q(x) 的最高项有贡献。
具体来说,递归到区间 [L,R] 的 Q(x) 只需要保留最高的 R-L+1 项。
这样时间复杂度就对了。
然后我们求出 Q(x) 后根据 R_i(x)=F(x)-Q_i(x)G_i(x) 算出 R(x) 即可。
基于转置原理的多项式多点求值算法
写得有点乱。
先令多项式项数和多点求值点数相同,均为 n ,若不同补齐即可。
对于一个 n-1 次多项式 \sum_{i=0}^{m-1}f_ix^i ,我们将其看做一个列向量 f=(f_0,\cdots,f_{m-1})^{T} 。那么对于 a_0,\cdots,a_{n-1} 多点求值的过程就相当于是计算:
\begin{bmatrix}
1&a_0&a_0^2&\cdots&a_0^{n-1}\\
1&a_1&a_1^2&\cdots&a_1^{n-1}\\
\vdots&\vdots&\vdots&&\vdots\\
1&a_{n-1}&a_{n-1}^2&\cdots&a_{n-1}^{n-1}\\
\end{bmatrix}
\cdot f
记左边的矩阵为 \mathbf{V} 。
考虑 \mathbf{V}^T\cdot f=g ,我们写出 g_i 的生成函数:
\sum_{i=0}^{n-1}g_ix^i=\sum_{i=0}^{n-1}x^i\sum_{j=0}^{n-1}a_j^if_j=\sum_{j=0}^{n-1}f_{j}\sum_{i=0}^{n-1}a_j^ix^i=\sum_{j=0}^{n-1}\frac{f_j}{1-a_jx}\bmod x^n
可以使用分治 NTT 计算,我们先求出和式的分子,再乘上和式的分母的逆即可。
但这对我们计算 \mathbf{V}\cdot f 有什么帮助呢?
这就要用到转置原理。其核心思想是:对于矩阵 \mathbf{V} 和任一向量 f ,为了优化计算 \mathbf{V}\cdot f ,我们可以先考察计算 \mathbf{V}^T\cdot f 的方法。然后我们把计算 \mathbf{V}^T\cdot f 的过程展开,看成是 \mathbf{E_k}\cdots \mathbf{E}_1\cdot f ,其中 \mathbf{E}_i 均为初等矩阵。那么 \mathbf{V}\cdot f s 的效果是使 f 的某个位置乘上 c ,那么 \mathbf{E}^T\cdot f 的效果也一样。
也就是说,我们把计算 \mathbf{V}^T\cdot f 的过程 “倒过来”,再对每一个操作 “转置” 一下,就可以得到计算 \mathbf{V} \cdot f 的操作过程了。
那么现在考虑分治 NTT 的过程,理论上由于它计算的是 \mathbf{V}^T\cdot f ,所以过程中应该是不断在对当前向量 u (初始为 f )做线性变换。可能会令人绕晕的一点是,分治 NTT 的过程中都是多项式在相乘,哪里出现了线性变换?但实质上是多项式的系数列表(即生成函数对应的原数列)在做线性变换:对于常多项式 A (次数为 m )和多项式 f (次数为 n )相乘,可以看做是在做类似于这样的变换:
\begin{bmatrix}A_0&0&\cdots\\A_1&A_0&0&\cdots\\\vdots&\vdots\\A_m&A_{m-1}&\cdots&A_0&0&\cdots\\0&A_m&\cdots&A_1&A_0&0&\cdots\\\vdots&\vdots\end{bmatrix}\begin{bmatrix}f_0\\f_1\\\vdots\\f_n\\0\\\vdots\end{bmatrix}=\begin{bmatrix}g_0\\g_1\\\vdots\\g_{n+m}\end{bmatrix}
那么我们会将 A_i f_j 贡献到 g_{i+j} ,即 g_{i+j}\mathrel{+}=A_i f_j ,注意这里应有 j\leq n ,因为在操作过程中我们是对整个向量的 f_0,\cdots,f_n 部分提取出来做多项式乘法,f_n 后面的那些位置并不是 0 。转置之后就是 f_j\mathrel{+}= A_ig_{i+j} ,其中需满足 j\leq n ,这恰好对应着减法卷积,记为 A \times^Tf 。(注意我们当然不可能真的把每次加减乘除交换操作全部丢到一个队列里面然后反过来做。我们要整体考虑一个操作的 “转置” 是什么)
考虑 NTT 中的过程:
分治求分子。假设当前分治到 [l,r] 。已经求得 [l,mid] 的答案 L=u_{l,mid} 和 [mid+1,r] 的答案 R=u_{mid+1,r} 。设 A_l=\prod_{i=l}^{mid}(1-a_ix),A_r=\prod_{i=mid+1}^r(1-a_ix) 。
那么当前的答案 F=A_r\times L+A_l\times R 。
已经求出了分子 P 。设 Q=\frac{1}{\prod_{i=0}^{n-1}(1-a_ix)}\bmod x^n ,执行操作 P\gets Q\times P \bmod x^n 。
那么整个过程转置之后就是:
执行操作 P\gets Q\times ^T P ,注意由于原来有取模所以这里转置后的多项式乘法会有所不同,具体表现在 f_j\mathrel{+}=A_i g_{i+j} 中还需满足 i+j<n (不过本来大于等于 n 也不会有项) 。然后将 P 扔进分治递归。
假设当前分治到 [l,r] 。原来我们进行的是 F=A_r\times L+A_l\times R (L=u_{l,mid},R=u_{mid+1,r} ),我们现在要把它转置回去。
可以看成是在同一个矩阵里同时进行了操作 F_{i+j}\mathrel{+}= A_{r,i}L_j 和 F_{i+j}\mathrel{+}=A_{l,i}R_j (注意两边的 j 都各有限制范围),那么转置回来就是 L_j\mathrel{+}=A_{r,i}F_{i+j} 和 R_j\mathrel{+}=A_{l,i}F_{i+j} 。即 L=A_r\times^T F,R=A_l\times^TF 。
然后递归下去做即可。
时间复杂度 O(n\log ^2n) ,常数较小。
多项式快速插值
洛必达法则
众所周知,两个无穷小之比或两个无穷大之比的极限可能存在,也可能不存在。因此,求这类极限时往往需要适当的变形,转化成可利用极限运算法则或重要极限的形式进行计算。洛必达法则便是应用于这类极限计算的通用方法。
——摘自百度百科
洛必达法则分为两种 0/0 型和 \infty/\infty 型:
若函数 $f(x),g(x)$ 满足以下条件:
- $\lim\limits_{x\to a}f(x)=0$,$\lim\limits_{x\to a}g(x)=0$。
- 在点 $a$ 的某去心邻域内两者都可导,且 $g'(x)\neq0$。
注意,这里只需要存在去心邻域即可,我们无需关心去心邻域的大小,因为我们只需要存在 $f'(x),g'(x)$($x\to a$),下面会讲。
则:
$$
\lim_{x\to a}\dfrac{f(x)}{g(x)}=\lim_{x\to a}\dfrac{f'(x)}{g'(x)}
$$
若函数 $f(x),g(x)$ 满足以下条件:
- $\lim\limits_{x\to a}f(x)=\infty$,$\lim\limits_{x\to a}g(x)=\infty$。
- 在点 $a$ 的某去心邻域内两者都可导,且 $g'(x)\neq0$。注意,这里同样只需要存在去心邻域即可,原因同上。
则:
$$
\lim_{x\to a}\dfrac{f(x)}{g(x)}=\lim_{x\to a}\dfrac{f'(x)}{g'(x)}
$$
证明:
考虑构造 u(x)=(g(x),f(x)) ,这个函数是一个 R\to R^2 的映射,即一条数轴映射到一个平面坐标系。
我们考虑求 u'(x) 是什么:
\begin{aligned}
u'(x)&=\lim_{\Delta x\to 0}\dfrac{f(x+\Delta x)-f(x)}{g(x+\Delta x)-g(x)}\\
&=\lim_{\Delta x\to 0}\dfrac{\quad\dfrac{f(x+\Delta x)-f(x)}{\Delta x}\quad}{\quad\dfrac{g(x+\Delta x)-g(x)}{\Delta x}\quad}\\
&=\lim_{\Delta x\to 0}\dfrac{f'(x)}{g'(x)}\\
&=\dfrac{f'(x)}{g'(x)}
\end{aligned}
对于一个非(0,0) 的点 A(x,y) ,我们考虑定义 A 的 “零点线” 为过 (0,0) 和 A 的直线,易知这条直线的斜率为 \dfrac{y}{x} 。
那么 u(x) 的零点线的斜率就是 \dfrac{f(x)}{g(x)} 。
显然我们现在已经构造出 \dfrac{f(x)}{g(x)} 和 \dfrac{f'(x)}{g'(x)} 的几何意义了,接下来就看什么时候他们相等了。
首先 \dfrac{f(x)}{g(x)} 为 u(x) 处的切线的斜率,\dfrac{f'(x)}{g'(x)} 为 u(x) 的零点线的斜率,如果他们相等,要么就是切线和零点线重合,要么就是切线和零点线平行。
切线和零点线重合(0/0 型):
大前提是 u(x) 过 (0,0) ,即存在 \lim\limits_{x\to a}f(x)=0 、\lim\limits_{x\to a}g(x)=0 的情况。
首先注意到切线和零点线都是过 u(x) 的,而零点线还过 (0,0) ,所以如果切线也过 (0,0) 那就能满足切线和零点线重合了。
易知当 u(x)\to(0,0) 时,切线肯定就是零点线。
所以有 \lim\limits_{u(x)\to(0,0)}\dfrac{f(x)}{g(x)}=\lim\limits_{u(x)\to(0,0)}\dfrac{f'(x)}{g'(x)} ,即为 0/0 型的洛必达法则。
切线和零点线平行(\infty/\infty 型):
大前提是 u(x) 过 (\infty,\infty) ,即存在 \lim\limits_{x\to a}f(x)=\infty 、\lim\limits_{x\to a}g(x)=\infty 的情况。
你可能会产生疑惑:既然切线和零点线都过 u(x) ,那么他们怎么可能平行呢?
但我们有这么一种特殊情况:相交于无穷远点的两条直线是平行的。
得到当 u(x)\to(\infty,\infty) 时,切线和零点线是趋于平行的,因为他们相交于一无穷远点 u(x) 。
但这个 “相交于无穷远点” 的解释对于我来说比较难以理解。
所以有 \lim\limits_{u(x)\to(\infty,\infty)}\dfrac{f(x)}{g(x)}=\lim\limits_{u(x)\to(\infty,\infty)}\dfrac{f'(x)}{g'(x)} ,即为 \infty/\infty 型的洛必达法则。
证毕。
(上述的证明比较感性,目的只是帮助你记忆,真正正确严格的证明请自行查阅资料)
多项式快速插值
先考虑用拉格朗日插值暴力搞出这个多项式:
f(x)=\sum_{i=1}^ny_i\prod_{j\neq i}\dfrac{x-x_j}{x_i-x_j}
显然这样做时间是无法接受的,考虑优化。
我们将它分为两个部分去计算:(因为分母的部分是可以直接计算出来的,接下来会讲)
f(x)=\sum_{i=1}^n\dfrac{y_i}{\prod\limits_{j\neq i}(x_i-x_j)}\prod_{j\neq i}(x-x_j)
先看到左边的那个式子,y_i 是个常数,我们暂时不用管它,接下来看分母:\prod\limits_{j\neq i}(x_i-x_j) 。
如果我们设 g(x)=\prod\limits_{i=1}^n(x-x_i) ,那么分母即为把 x=x_i 代入 \dfrac{g(x)}{x-x_i} 得到的值。但代入后分子分母都是 0 ,怎么求啊……
这就要用到 0/0 型洛必达法则了:
\lim_{x\to x_i}\dfrac{g(x)}{x-x_i}=\lim_{x\to x_i}\dfrac{g'(x)}{(x-x_i)'}=\lim_{x\to x_i}g'(x)=g'(x_i)
那么我们只需要用分治 NTT 求出 g(x) ,再求导得出 g'(x) ,再代入 x_1,x_2,\cdots,x_n 多点求值即可。
考虑分治计算最后的式子,定义 f_{l,r}(x) 为:
\begin{aligned}
f_{l,r}(x)&=\sum_{i=l}^r\dfrac{y_i}{\prod\limits_{j=1,j\neq i}^n(x_i-x_j)}\prod_{j=l,j\neq i}^r(x-x_j)\\
&=\sum_{i=l}^r\dfrac{y_i}{g'(x_i)}\prod_{j=l,j\neq i}^r(x-x_j)
\end{aligned}
容易得到:
f_{l,r}(x)=f_{l,mid}(x)\left[\prod_{j=mid+1}^r(x-x_j)\right]+f_{mid+1,r}(x)\left[\prod_{j=l}^{mid}(x-x_j)\right]
然后就可以分治求了。
拉格朗日反演
对于形式幂级数(以下简称幂级数) F(x) ,记 F[i]=[x^i]F(x) 。
幂级数复合:(F\circ G)(x)=F(G(x))=\sum_{i\geq 0}F[i]G(x)^i 。可以看成一个幂级数变换。
常数项为 0 ,一次项非 0 的幂级数集合 \mathbf{P} 在复合运算下形成群:
封闭性:对于任意 F,G\in \mathbf{P} ,(F\circ G)(x) 常数项为 0 ,一次项系数为 F[1]G[1]\neq 0 。
结合律:映射的复合显然有结合律。
单位元:F(x)=x 。
逆元:即对于任意的 F\in \mathbf{P} ,都存在 F^{<-1>}\in \mathbf{P} 满足 (F^{<-1>}\circ F)(x)=x 。根据幂级数复合的定义 (F^{<-1>}\circ F)(x)=\sum_{i\geq 0}F^{<-1>}[i]F(x)^i=x 递推构造即可得到唯一的 F^{<-1>} 。
补充的几点:
注意幂级数复合并不满足交换律,但满足 (F^{<-1>}\circ F)(x)=(F\circ F^{<-1>})(x)=x ,称 F,F^{<-1>} 互为复合逆。
为什么一定要求常数项为 0 、一次项非 0 才能形成群?
接下来,我们引入分式域,即引入 x 的负整数次幂这种东西。接下来我们对于幂级数的操作都在分式域上进行。
分式域上的幂级数加减乘仍然是良定义的。而对于分式域上的逆,考虑要求 F(x) 的逆,我们将 F(x) 进行适当的平移使得 F(x) 变为整式,即找到一个整数 k 使得 F(x)=x^kG(x) 且 G(x) 为常数项非 0 的整式,就能得到 F(x)^{-1}=x^{-k}G(x)^{-1} 。
分式域上的求导、积分仍然是良定义的,不过有一个等会要用到的性质就是:对于任意分式域上的幂级数 F(x) ,[x^{-1}]F(x)'=0 。
引理 1 :对于任意的 F(x)\in \mathbf{P} ,有:
[x^{-1}]F'(x)F(x)^k=[k=-1]
证明 :当 k\neq -1 时,F'(x)F(x)^k=(\frac{1}{k+1}F(x)^{k+1})' ,再根据刚刚提到的性质即可得证。
当 k=-1 时,[x^{-1}]F'(x)F(x)^k=[x^{-1}]\frac{F'(x)}{F(x)}=[x^0]\frac{F'(x)}{F(x)/x}=\frac{F[1]}{F[1]}=1 。
拉格朗日反演 :对于 F(x),G(x)\in\mathbf{P} ,若 F,G 互为复合逆,那么有:
[x^n]G(x)=\frac{1}{n}[x^{-1}]F(x)^{-n}
证明 :考虑利用引理1把 G[n] 提出来:
\begin{aligned}
G(F)&=x\\
\sum_{i\geq 0}G[i]F^{i}&=x\\
\sum_{i\geq 0}G[i]F^{i-n-1}F'&=xF^{-n-1}F'\\
[x^{-1}]\sum_{i\geq 0}G[i]F^{i-n-1}F'&=[x^{-1}]xF^{-n-1}F'\\
\sum_{i\geq 0}G[i]\cdot [i-n-1=-1]&=[x^{-1}]xF^{-n-1}F'\\
G[n]&=[x^{-1}]xF^{-n-1}F'\\
G[n]&=[x^{-2}]-\frac{1}{n}(F^{-n})'\\
G[n]&=\frac{1}{n}[x^{-1}]F^{-n}
\end{aligned}
计算方法 :一个更方便的形式是 \frac{1}{n}[x^{n-1}](F(x)/x)^{-n} ,于是我们可以直接先求逆然后做多项式快速幂就能 O(n\log n) 计算 G[n] 了。
扩展拉格朗日反演 :对于 F(x)\in\mathbf{P} 和另外一个无要求的幂级数 H ,若 T(F(x))=H(x) ,F,G 互为复合逆,有:
[x^n]T(x)=[x^n]H(G(x))=\frac{1}{n}[x^{-1}]H'(x)F(x)^{-n}
证明 :也是类似的:
\begin{aligned}
T(F)&=H\\
\sum_{i\geq 0}T[i]F^{i}&=H\\
[x^{-1}]\sum_{i\geq 0}T[i]F^{i-n-1}F'&=[x^{-1}]H\cdot F^{-n-1}F'\\
T[n]&=[x^{-1}]H\cdot F^{-n-1}F'\\
T[n]&=[x^{-1}]-\frac{1}{n}H\cdot (F^{-n})'\\
T[n]&=-\frac{1}{n}\sum_i H[i]\times (F^{-n})'[-1-i]\\
T[n]&=-\frac{1}{n}\sum_i(-i) H[i]\times F^{-n}[-i]\\
T[n]&=\frac{1}{n}\sum_i H'[i-1]\times F^{-n}[-i]\\
T[n]&=\frac{1}{n}[x^{-1}]H'\times F^{-n}\\
\end{aligned}
计算方法 :也是类似的变成 [x^{n-1}]H'\times(F(x)/x)^{-n} 即可 O(n\log n) 求解单项 T[n] 。
应用 :目前遇到了两种:
例:要求计算 [x^n]B(x)^k ,其中 B(x)=\dfrac{1-\sqrt{1-4x}}{2x} 为卡特兰数生成函数。
这里显然是设 H(x)=x^k ,但我们不能直接对 B(x) 求符合逆,因为它常数项不为 0 ,所以我们令 B'(x)=\frac{1-\sqrt{1-4x}}{2} ,变为计算 [x^n]x^{-k}H(B'(x))=[x^{n+k}]H(B'(x)) ,其中 B'(x)\in \mathbf{P} 。
接下来是计算 B'(x) 的复合逆 A(x) :
\begin{aligned}
y&=\frac{1-\sqrt{1-4x}}{2}\\
1-2y&=\sqrt{1-4x}\\
4y^2-4y+1&=1-4x\\
x&=y-y^2
\end{aligned}
于是得到 A(x)=x-x^2 ,根据扩展拉格朗日反演:
\begin{aligned}
\ [x^{n}]H(B'(x)))&=\frac{1}{n}[x^{n-1}]H'(x)\left(\frac{x}{A(x)}\right)^{n}\\
&=\frac{1}{n}[x^{n-1}]kx^{k-1}\left(\frac{1}{1-x}\right)^{n}\\
&=\frac{k}{n}[x^{n-k}]\left(\frac{1}{1-x}\right)^{n}\\
&=\frac{k}{n}\binom{2n-k-1}{n-1}\\
[x^{n+k}]H(B'(x))&=\frac{k}{n+k}\binom{2n+k-1}{n+k-1}
\end{aligned}
FWT
位运算卷积(FWT)
子集卷积
对于数组 A,B ,求数组 C 使得 C_k=\sum\limits_{i|j=k,i\&j=0}A_iB_j 。
直观地看上去就是 C_k=\sum\limits_{i\in k}A_iB_{k-i} ,所以称作子集卷积。
下面用 |i| 表示 \operatorname{popcount}(i) 。
注意到 i|j=k,i\&j=0 等价于 i|j=k,|i|+|j|=|k| 。
我们设:
a_{i,j}=
\begin{cases}
A_j&\text{if }i=|j|\\
0&\text{if }i\neq|j|
\end{cases}
对于 b_{i,j} 的定义同理。
那我们要求的实际是:
C_k=\sum_{i|j=k,|i|+|j|=|k|}A_iB_j\\
=\sum_{p+q=|k|}\sum_{i|j=k}a_{p,i}b_{q,j}\\
=\sum_{p+q=|k|}\operatorname{OR}(a_p,b_q)_k\\
=\left[\sum_{p+q=|k|}\operatorname{OR}(a_p,b_q)\right]_k
于是我们设 c_{r,k}=\left[\sum_{p+q=r}\operatorname{OR}(a_p,b_q)\right]_k ,那么我们最后要求的答案是 C=\{c_{|i|,i}\} 。
那么我们可以先 O(\log n) 枚举 p ,把每个 a_p 的 FWT 值给算出来。同理我们也把每个 b_q 的 FWT 的值都算出来。然后再 O(\log^2 n) 枚举 p,q ,把每个 c_r 的 FWT 值给算出来。最后再 O(\log n) 枚举 r ,做一遍 IFWT 把 c_r 给逆推回来,那么也就得到了 C 。
总共的时间复杂度为 O(n\log^2 n) 。
个人认为这个算法的思想是把 O(n) 的条件 i\&j=k 转化成 O(\log n) 的条件 |i|+|j|=|k| ,最后再利用 FWT 把时间复杂度从 O(n^2\log^2 n) 降到了 O(n\log^2n) 。
代码可参考【XSY3020】CF914G Sum the Fibonacci。
一些特殊形式的求法
一、求 \sum\limits_{k=1}^{n}x^k\sum\limits_{i-j=k}a_ib_j
平常的 FFT 是求 \sum\limits_{k=1}^{n}x^k\sum\limits_{i+j=k}a_ib_j ,那么我现在让你求 \sum\limits_{k=1}^{n}x^k\sum\limits_{i-j=k}a_ib_j 。
设 S(k)=\sum\limits_{i-j=k}a_ib_j ,那么题目是要求 \sum\limits_{k=1}^nx^kS(k) 。
设 \widehat{a_i}=a_{n-i} ,那么:
\begin{aligned}
S(k)=&\sum\limits_{i-j=k}a_ib_j\\
=&\sum_{(n-i)-j=k}\widehat{a_i}b_j\\
=&\sum_{i+j=n-k}\widehat{a_i}b_j
\end{aligned}
那么:
S(n-k)=\sum_{i+j=n-(n-k)}\widehat{a_i}b_j=\sum_{i+j=k}\widehat{a_i}b_j
那么就可以设多项式 A=\{\widehat{a_i}\} ,B=\{b_i\} ,C=A\times B ,则 S(n-k)=[x^k]C 。
所以将 C 翻转一下就是我们要求的了。
拉格朗日插值
原理
给出 n 个点 (x_i,y_i) ,可知这 n 个点唯一确定一个多项式 f(x) 。
现在询问 k ,求 f(k) 。
用高斯消元可以做到预处理 O(n^3) ,每次询问 O(n) 。
但是现在只有一次询问,且 n 很大。
所以我们考虑用询问更快的方法求解。
拉格朗日插值是通过构造的方法,在 O(n^2) 的时间内进行询问。
设拉格朗日基本多项式为:
\ell_i(x)=\prod_{j=1,j\neq i}^{n}\frac{x-x_j}{x_i-x_j}
注意到当 x=x_i 时,\ell_i(x)=1 ;当 x=x_j (j\neq i )时,\ell_i(x)=0 。
那么就能把这个多项式 f(x) 构造出来了:
f(x)=\sum_{i=1}^n y_i\ell_i(x)
显然,这个多项式满足对于任意的 i ,f(x_i)=y_i 。
那么这个多项式就是我们要找的多项式。
\min-\max 容斥
\min-\max 容斥
设 S 为一集合,\min(S) 表示集合 S 的最小元素,\max(S) 表示集合 S 的最大元素。有:
\max(S)=\sum_{T\subseteq S,T\neq \emptyset}(-1)^{|T|-1}\min(T)
**$\text{kth}\min-\max$ 容斥**
设 $S$ 为一集合,$\min_k(S)$ 表示集合 $S$ 的第 $k$ 小元素,$\max_k(S)$ 表示集合 $S$ 的第 $k$ 大元素。有:
$$
\max\nolimits_k(S)=\sum_{T\subseteq S,T\neq \emptyset}(-1)^{|T|-k}\binom{|T|-1}{k-1}\min(T)
$$
改成 $\min_k$ 和 $\max$ 相关的式子也是一样的。
证明:
下面给出的是 $\text{kth}\min-\max$ 容斥的证明,而 $\min-\max$ 容斥的证明类似。
先设容斥系数 $f(x)$ 使得:
$$
\max\nolimits_k(S)=\sum_{T\subseteq S,T\neq \emptyset}f(|T|)\min(T)
$$
设 $n=|S|$,考虑 $S$ 的第 $x$ 小的数 $v$ 在式子左右两边的贡献:
- 在左边贡献为 $v[x=n-k+1]$。
- 在右边,我们枚举所有包含 $v$ 且以 $v$ 作为最小值的 $T$,贡献为 $\sum\limits_{i=0}^{n-x}\dbinom{n-x}{i}f(i+1)v$。
于是:
$$
\begin{aligned}
v[x=n-k+1]&=\sum\limits_{i=0}^{n-x}\dbinom{n-x}{i}f(i+1)v\\
[x=n-k+1]&=\sum\limits_{i=0}^{n-x}\dbinom{n-x}{i}f(i+1)\\
[n-x=k-1]&=\sum\limits_{i=0}^{n-x}\dbinom{n-x}{i}f(i+1)\\
\end{aligned}
$$
设 $F(x)=f(x+1)$,$G(x)=[x=k-1]$,令 $n'=n-x$,于是:
$$
\begin{aligned}
G(n')&=\sum\limits_{i=0}^{n'}\dbinom{n'}{i}F(i)\\
\end{aligned}
$$
根据二项式反演:
$$
\begin{aligned}
F(n)&=\sum\limits_{i=0}^{n}(-1)^{n-i}\dbinom{n}{i}G(i)\\
\end{aligned}
$$
带回去:
$$
\begin{aligned}
f(n+1)&=\sum_{i=0}^n(-1)^{n-i}\binom{n}{i}[i=k-1]\\
&=(-1)^{n-k+1}\binom{n}{k-1}
\end{aligned}
$$
故:
$$
f(n)=(-1)^{n-k}\binom{n-1}{k-1}
$$
证毕。
上面最开始的两个式子套上期望也是成立的。口胡证明:因为期望=∑(枚举每种情况)概率×贡献,而对于每种情况都适用于上面不带期望的式子,再根据期望的线性性可推得成立。
$$
\mathbb{E}(\max(S))=\sum_{T\subseteq S,T\neq \emptyset}(-1)^{|T|-1}\mathbb{E}(\min(T))
$$
$$
\mathbb{E}(\max\nolimits_k(S))=\sum_{T\subseteq S,T\neq \emptyset}(-1)^{|T|-k}\binom{|T|-1}{k-1}\mathbb{E}(\min(T))
$$
## 斯特林数
### 第二类斯特林数
第二类斯特林数 $\begin{Bmatrix}n\\m\end{Bmatrix}$ 表示将 $n$ 个不同的球分到 $m$ 个没有区别的袋子里,每个袋子里不能为空,且每个袋子里的球无序的方案数。
#### 递推公式
$$
\begin{aligned}
\begin{Bmatrix}n\\0\end{Bmatrix}&=[n=0]\\
\begin{Bmatrix}n\\m\end{Bmatrix}&=\begin{Bmatrix}n-1\\m-1\end{Bmatrix}+m\begin{Bmatrix}n-
1\\m\end{Bmatrix}
\end{aligned}
$$
证明:从组合意义证明,我们考虑前 $n-1$ 个球已经装好了,对于最后一个球来说,它有可能单独新放进一个空的袋子,也可能在非空的袋子中任选一个。
#### 通项公式
$$
\begin{Bmatrix}n\\m\end{Bmatrix}=\dfrac{1}{m!}\sum_{i=0}^m(-1)^i\dbinom{m}{i}(m-i)^n
$$
感性理解是:容斥,先枚举有 $i$ 个袋子是空的,再枚举是哪 $i$ 个袋子,然后 $n$ 个球就随便放进剩下的 $m-i$ 的袋子即可。又由于袋子是无序的,所以需要除以 $m!$。
拆开组合数,再往下推下去:
$$
\begin{aligned}
\begin{Bmatrix}n\\m\end{Bmatrix}&=\dfrac{1}{m!}\sum_{i=0}^m(-1)^i\dbinom{m}{i}(m-i)^n\\
&=\sum_{i=0}^m\dfrac{(-1)^i(m-i)^n}{i!(m-i)!}
\end{aligned}
$$
#### 性质
$$
m^n=\sum_{i=0}^{\min(n,m)}\binom{m}{i}\begin{Bmatrix}n\\i\end{Bmatrix}i!
$$
从组合意义解释:
- 等式左边的 $m^n$ 表示将 $n$ 个不同的小球扔到 $m$ 个不同的盒子里,盒子可以为空的方案数。
- 等式右边先枚举非空的盒子个数 $i$(注意非空盒子个数不能大过总盒子数,也不能大过总球数,故 $i$ 的枚举上界为 $\min(n,m)$),再从这 $m$ 个盒子中选出这 $i$ 个盒子,然后把 $n$ 个球扔到这 $i$ 个盒子里,不能留空,最后再乘上 $i!$,因为 $\begin{Bmatrix}n\\i\end{Bmatrix}$ 的定义中盒子是相同的,但现在盒子是不同的。
事实上,$i$ 的枚举上界是 $n$ 是 $m$ 都可以,因为当 $i>\min(n,m)$ 时,$\dbinom{m}{i}$ 和 $\begin{Bmatrix}n\\i\end{Bmatrix}$ 必有一个为 $0$。
再往下推下去还能跟下降幂扯上关系:
$$
m^n=\sum_{i=0}^{\min(n,m)}\binom{m}{i}\begin{Bmatrix}n\\i\end{Bmatrix}i!=\sum_{i=0}^n\begin{Bmatrix}n\\i\end{Bmatrix}m^{\underline{i}}
$$
#### 第二类斯特林数·行
给定 $n$,对于所有的 $i\in [0,n]$,求出 $\begin{Bmatrix}n\\i\end{Bmatrix}$。
注意到上面通项公式有 $\begin{Bmatrix}n\\m\end{Bmatrix}=\sum\limits_{i=0}^m\dfrac{(-1)^i(m-i)^n}{i!(m-i)!}$,那么我们直接设 $A_i=\dfrac{(-1)^i}{i!}$、$B_i=\dfrac{i^n}{i!}$,然后多项式相乘 $A\times B$ 即得答案。
#### 第二类斯特林数·列
给定 $n,k$,对于所有的 $i\in [0,n]$,求出 $\begin{Bmatrix}i\\k\end{Bmatrix}$。
根据第二类斯特林数的递推式:
$$
\begin{Bmatrix}n\\m\end{Bmatrix}=\begin{Bmatrix}n-1\\m-1\end{Bmatrix}+m\begin{Bmatrix}n-
1\\m\end{Bmatrix}
$$
设 $F_k(x)=\sum\limits_{i=0}\begin{Bmatrix}i\\k\end{Bmatrix}x^i$,那么:
$$
\begin{aligned}
F_k(x)&=\sum\limits_{i=0}\begin{Bmatrix}i\\k\end{Bmatrix}x^i=\sum\limits_{i=1}\begin{Bmatrix}i\\k\end{Bmatrix}x^i\\
&=\sum_{i=1}\left(\begin{Bmatrix}i-1\\k-1\end{Bmatrix}+k\begin{Bmatrix}i-1\\k\end{Bmatrix}\right)x^i\\
&=\sum_{i=1}\begin{Bmatrix}i-1\\k-1\end{Bmatrix}x^i+\sum_{i=1}k\begin{Bmatrix}i-1\\k\end{Bmatrix}x^i\\
&=\sum_{i=1}\begin{Bmatrix}i\\k-1\end{Bmatrix}x^{i+1}+\sum_{i=0}k\begin{Bmatrix}i\\k\end{Bmatrix}x^{i+1}\\
&=xF_{k-1}(x)+kxF_k(x)
\end{aligned}
$$
解得 $F_k(x)=\dfrac{x}{1-kx}F_{k-1}(x)$。
容易得到边界 $F_0(x)=1$,那么:
$$
F_k(x)=\dfrac{x}{1-kx}F_{k-1}(x)=\cdots=\dfrac{x^k}{\prod\limits_{i=1}^k(1-ix)}
$$
分治 NTT 求分母,然后再求逆即可。时间复杂度 $O(n\log ^2n)$。
当然还有时间复杂度为 $O(n\log n)$ 的更优秀做法,实现方法和第一类斯特林数·列类似。
### 第一类斯特林数
第一类斯特林数 $\begin{bmatrix}n\\m\end{bmatrix}$ 表示将 $n$ 个不同的球划分为 $k$ 个互不区分且不能为空的轮换的方案数。
轮换可以理解为环,旋转后相同的轮换互不区分,顺时针和逆时针的轮换需要区分(如 $[A,B,C,D]\neq [D,C,B,A]$)。
#### 递推公式
$$
\begin{aligned}
\begin{bmatrix}n\\0\end{bmatrix}&=|n=0|
\\\begin{bmatrix}n\\m\end{bmatrix}&=\begin{bmatrix}n-1\\m-1\end{bmatrix}+(n-1)\begin{bmatrix}n-1\\m\end{bmatrix}
\end{aligned}
$$
证明:从组合意义证明,我们考虑前 $n-1$ 个球已经分好了,现在插入第 $n$ 个球。若前面的球组成了 $m-1$ 个轮换,那么这个球就得自己组成一个轮换;若前面的球已经组成了 $m$ 个轮换,那么这个球可以跟在前 $n-1$ 个球的任意一个球后面。
第一类斯特林数没有实用的通项公式。
#### 第一类斯特林数·行
结论:$\begin{bmatrix}n\\i\end{bmatrix}$ 的一般生成函数(OGF)为 $x^{\overline{n}}$。
证明:使用数学归纳法。
$$
\begin{aligned}
x^{\overline{n+1}}&=(x+n)x^{\overline{n}}\\
&=(x+n)\sum_{i=0}^n\begin{bmatrix}n\\i\end{bmatrix}x^i\\
&=\sum_{i=0}^n\left(\begin{bmatrix}n\\i-1\end{bmatrix}+n\begin{bmatrix}n\\i\end{bmatrix}\right)x^i\\
&=\sum_{i=0}^n\begin{bmatrix}n+1\\i\end{bmatrix}x^i\\
\end{aligned}
$$
于是我们求出 $x^{\overline{n}}$ 即可,$x^{\overline{n}}=\prod\limits_{i=0}^{n-1}(x+i)$。
直接分治 NTT 是 $O(n\log ^2n)$ 的,不够优秀。
考虑倍增:假设我们已经求出了 $x^{\overline{n}}$,现在要求 $x^{\overline{2n}}$。
设 $f(x)=x^{\overline{n}}$,那么 $x^{\overline{2n}}=f(x)f(x+n)$。现在问题变为了:已经知道了 $f(x)$,如何求 $f(x+n)$?
我们直接暴力展开就好:
$$
\begin{aligned}
f(x+n)&=\sum_{i=0}^nf_i(x+n)^i\\
&=\sum_{i=0}^nf_i\sum_{j=0}^i\binom{i}{j}n^{i-j}x^j\\
&=\sum_{i=0}^nf_i\sum_{j=0}^i\dfrac{i!}{j!(i-j)!}n^{i-j}x^j\\
&=\sum_{j=0}^n\dfrac{1}{j!}x^j\sum_{i=j}^{n}f_ii!\times \dfrac{n^{i-j}}{(i-j)!}\\
\end{aligned}
$$
那么就可以卷积做了。
于是从 $x^{\overline{n}}$ 推到 $x^{\overline{2n}}$ 需要先 $O(n\log n)$ 求出 $f(x+n)$,再把 $f(x)$ 和 $f(x+n)$ 相乘 $O(n\log n)$ 得到 $x^{\overline{2n}}$。
于是总时间复杂度为 $O(n\log n)$。
#### 第一类斯特林数·列
待填坑……
## 其他
### 自然数幂求和
求 $S_k(n)=\sum\limits_{i=0}^ni^k$。
#### 1. $n$ 固定,询问 $k=0\sim K$ 时的答案。
**算法特点:时间复杂度 $O(K\log K)$,需要模数能支持 NTT。**
我们求出一个与 $S_k(n)$ 有关的生成函数:
$$
\begin{aligned}
&\sum_{k=0}^KS_k(n)\frac{1}{k!}x^k\\
=&\sum_{k=0}^K\sum_{i=0}^n i^k\frac{1}{k!}x^k\\
=&\sum_{i=0}^n\sum_{k=0}^K \frac{(ix)^k}{k!}
\end{aligned}
$$
注意到 $e^x$ 的泰勒展开为:$e^x=1+\dfrac{x}{1!}+\dfrac{x^2}{2!}+\cdots$,故原式为:
$$
=\sum_{i=0}^ne^{ix}\pmod {x^{K+1}}
$$
注意到 $\dfrac{1}{1-x}=1+x+x^2+\cdots$,有 $1+x+x^2+\cdots+x^n=\dfrac{1-x^{n+1}}{1-x}$。
那么原式即为:
$$
\begin{aligned}
=&\sum_{i=0}^n(e^x)^i\pmod {x^{K+1}}\\
=&\frac{1-e^{(n+1)x}}{1-e^x}\pmod {x^{K+1}}\\
=&\frac{e^{(n+1)x}-1}{e^x-1}\pmod {x^{K+1}}\\
=&\dfrac{\dfrac{(n+1)x}{1!}+\dfrac{(n+1)^2x^2}{2!}+\dfrac{(n+1)^3x^3}{3!}+\cdots}{\dfrac{x}{1!}+\dfrac{x^2}{2!}+\dfrac{x^3}{3!}+\cdots}\pmod {x^{K+1}}\\
=&\dfrac{(n+1)+\dfrac{(n+1)^2}{2!}x+\dfrac{(n+1)^3}{3!}x^2+\cdots}{1+\dfrac{1}{2!}x+\dfrac{1}{3!}x^2+\cdots}\pmod {x^{K+1}}\\
=&\dfrac{(n+1)+\dfrac{(n+1)^2}{2!}x+\cdots+\dfrac{(n+1)^{K+1}}{(K+1)!}x^{K}}{1+\dfrac{1}{2!}x+\cdots+\dfrac{1}{(K+1)!}x^{K}}\\
\end{aligned}
$$
多项式求逆即可。
最后得到的系数记得乘上个 $k!$ 才是 $S_k(n)$。
时间复杂度 $O(K\log K)$。
~~听说这玩意叫伯努利数?~~
#### 2. 对于特定的 $n,k$。
**算法特点:时间复杂度近似 $O(k)$,需要模数能支持求逆。**
首先证明 $S_k(n)$ 是关于 $n$ 的 $k+1$ 次多项式:
考虑序列 $a$,其中 $a_n=S_k(n)=\sum\limits_{i=1}^n i^k$。
那么序列 $a$ 的一阶差分序列 $b_n=a_{n+1}-a_n=\sum\limits_{i=1}^{n+1} i^k-\sum\limits_{i=1}^n i^k=(n+1)^k$,是一个关于 $n$ 的 $k$ 次多项式。
那么序列 $b$ 就是一个 $k$ 阶等差数列,那么序列 $a$ 的一阶差分序列是一个 $k$ 阶等差数列,那么序列 $a$ 是一个 $k+1$ 阶等差数列,所以 $a_n$ 的通项公式是一个关于 $n$ 的 $k+1$ 次多项式。
证毕。
但是我们不知道这个 $k+1$ 多项式到底长啥样,而又注意到我们只需要求代入 $n$ 时的多项式的值,所以考虑拉格朗日插值。
又有一个好处就是我们要代入的 $k+2$ 个点可以自己选。
我们考虑代入点集 $\{(1,S_k(1)),(2,S_k(2)),\cdots,(k+2,S_k(k+2))\}$。
那么我们就是要求出:
$$
\begin{aligned}
f(n)=&\sum_{i=1}^{k+2}S_k(i)\prod_{j\neq i}\frac{n-x_j}{x_i-x_j}\\
=&\sum_{i=1}^{k+2}S_k(i)\prod_{j\neq i}\frac{n-j}{i-j}\\
\end{aligned}
$$
先观察右边的 $\prod\limits_{j\neq i}\dfrac{n-j}{i-j}$:
- 对于分母,我们可以 $O(k)$ 预处理出 $pre_i=\prod\limits_{j=1}^i(n-j)$ 和 $suf_i=\prod\limits_{j=i}^n(n-j)$。
- 对于分子,我们可以 $O(k)$ 预处理出阶乘 $fac_i=i!$ 和它们的逆元 $inv_i=\dfrac{1}{fac_i}$。
那么 $\prod\limits_{j\neq i}\dfrac{n-j}{i-j}=\dfrac{pre_{i-1}suf_{i+1}}{(-1)^{k+2-i}fac_ifac_{k+2-i}}=(-1)^{k+2-i}pre_{i-1}\cdot suf_{i+1}\cdot inv_i\cdot inv_{k+2-i}$。
现在的问题就是如何对于每一个 $i\in[1,k+2]\cap\mathbb{Z}$ 求出 $S_k(i)$。也就是如何对于每一个 $i$ 求出 $i^k$,然后求 $S_k(i)$ 时只需求一遍前缀和即可。
如果直接暴力每次做一遍快速幂的话,时间复杂度为 $O(k\log k)$。
注意到 $i^k$ 为完全积性函数,所以我们直接线性筛即可。
那么此时只需对所有 $\leq k+2$ 的质数 $p$ 求 $p^k$。不妨设 $\leq k$ 的质数个数为 $\pi(k)$,那么时间复杂度为 $O(\pi(k)\log k)$。
又由于 $\pi(k)\approx \dfrac{k}{\ln(k)}$,$\ln(k)\approx \log(k)$,故时间复杂度近似为 $O(k)$。(实测当 $k=10^6$ 时 $\pi(k)\log k\approx 1.5\times 10^6$)。
#### 3. 使用第二类斯特林数
**算法特点:不需要模数支持取逆。**
根据第二类斯特林数的性质:
$$
m^n=\sum_{i=0}^{n}\binom{m}{i}\begin{Bmatrix}n\\i\end{Bmatrix}i!=\sum_{i=0}^n\begin{Bmatrix}n\\i\end{Bmatrix}m^{\underline{i}}
$$
于是:
$$
\begin{aligned}
&\sum_{i=0}^ni^k\\
=&\sum_{i=0}^n\sum_{j=0}^{k}\begin{Bmatrix}k\\j\end{Bmatrix}i^{\underline j}\\
=&\sum_{j=0}^k\begin{Bmatrix}k\\j\end{Bmatrix}\sum_{i=0}^ni^{\underline j}\\
=&\sum_{j=0}^k\begin{Bmatrix}k\\j\end{Bmatrix}\dfrac{(n+1)^{\underline{j+1}}}{j+1}
\end{aligned}
$$
最后一步证明可以用二项式反演:$\sum_{i=0}^ni^{\underline{j}}=j!\sum_{i=0}^n\binom{i}{j}=j!\binom{n+1}{j+1}=\frac{(n+1)^{\underline{j+1}}}{j+1}$。
注意这里看上去除了个 $j+1$,但我们不需要在模意义下求逆,因为 $n-j+1,n-j+1,\cdots,n+1$ 这 $j+1$ 个连续自然数肯定恰好有一个是 $j+1$ 的倍数。