我们知道一元二次方程有求根公式可以轻松的求根，一元三次方程也有，但稍微有些复杂，一元四次方程也有一个复杂的求根公式，但到了一元五次方程开始，除了特殊构造的方程外，大部分没有求根公式可以求出根，这时候就需要一些数值方法来逼近根，这些方法大多用于高于2次的方程中，这里我要介绍的就是一种名为 bairstow's method。不过我的方法稍微和原来的不太一样。
该方法思路很简单，首先将多项式：
$$ f(x) = sum_{i = 0}^{N} a_i x^i $$
分解为一个一元二次方程和一个 $ N - 2 $ 次多项式的乘积以及一个余项和常数项（注意这里有个地方不一样，我自己处理的时候要求 $ b_{N-2}x^{N-2} $ 这项的系数 $ b_{N-2} $一定为 $ 1 $，因为我会把原方程 $ f(x) $ 处理为 $ f(x) = frac{f(x)}{a_n} $，这样也会简单很多），即：
$$ f(x) = (x^2 + ux + v)(sum_{j=0}^{N-2} b_jx^j) + cx + d $$
如果我们将该方程还原，然后对照原多项式齐次项，把系数写出来，可以得到 $ c $、$ d $、$ a_i $、$ b_j $、$ u $、$ v $ 的关系，嘛，实际上就会得到：
$$ egin{aligned} c &= a_1 - b_0u - b_1v \ d &= a_0 - b_0v \ b_i &= a_{i+1} - b{i+1}u - b{i+2}v \ &... \ a_n &= b_{n-2} \ b_n &= b_{n - 1} = 0 end{aligned} $$
按照我的处理的话，这里 $ a_n = b_{n - 2} = 1 $，然后我们就只需要重点关注 $ c $、$ d $ 了。
考虑分解后的方程，我们令 $ f(x) = 0 $，这样我们就只需要考虑如何选取 $ u $、$ v $ 使得 $ c = 0 $、$ d = 0 $，考虑引入二维牛顿-拉弗森法
$$ egin{aligned} c(x, y) - c(u, v) &= frac{partial c}{partial u} (x - u) + frac{partial c}{partial v} (y - v) \ d(x, y) - d(u, v) &= frac{partial d}{partial u} (x - u) + frac{partial d}{partial v} (y - v) end{aligned} $$
令 $ c(x, y) = d(x, y) = 0 $，再用矩阵形式来重写式子：
$$ -egin{bmatrix} c \ d end{bmatrix} = egin{bmatrix}frac{partial c}{partial u} & frac{partial c}{partial v} \ frac{partial d}{partial u} & frac{partial d}{partial v} end{bmatrix} cdot egin{bmatrix} x \ y end{bmatrix} - egin{bmatrix} frac{partial c}{partial u} & frac{partial c}{partial v}\ frac{partial d}{partial u} & frac{partial d}{partial v} end{bmatrix} cdot egin{bmatrix} u \ v end{bmatrix} $$
然后等式两边同时乘以一个 $ egin{bmatrix} frac{partial c}{partial u} & frac{partial c}{partial v} \ frac{partial d}{partial u} & frac{partial d}{partial v} end{bmatrix}^{-1} $ ，得到：
$$ egin{aligned} -egin{bmatrix} frac{partial c}{partial u} & frac{partial c}{partial v}\ frac{partial d}{partial u} & frac{partial d}{partial v}end{bmatrix}^{-1}egin{bmatrix}c \ dend{bmatrix} &= egin{bmatrix}x \ yend{bmatrix} - egin{bmatrix}u \ vend{bmatrix} \egin{bmatrix} u \ v end{bmatrix}-egin{bmatrix} frac{partial c}{partial u} & frac{partial c}{partial v}\ frac{partial d}{partial u} & frac{partial d}{partial v} end{bmatrix}^{-1}egin{bmatrix}c \ d end{bmatrix} &= egin{bmatrix} x \ y end{bmatrix} end{aligned} $$
再令 $ u = u_i $、$ v = v_i $、 $ x = u_{i+1} $、$ y = v_{i+1} $，得到递推式：
$$ egin{bmatrix} u_{i+1} \ v_{i+1}end{bmatrix} = egin{bmatrix} u_i \ v_i end{bmatrix}-egin{bmatrix} frac{partial c}{partial u} & frac{partial c}{partial v}\ frac{partial d}{partial u} & frac{partial d}{partial v} end{bmatrix}^{-1} egin{bmatrix} c \ d end{bmatrix} $$
于是，当 $ c $、$ d $ 无限趋近于0的时候，原式变为 $ (x^2+ux + v)(sum_{i=0}^{N-2}b_ix^i) = 0 $，这样就可以使用一元二次方程求根公式得到等式左边第一部分，如果 $ N-2 > 2 $ 我们将第二部分递归使用该算法即可将所有根找到，当然这里还有一个重要部分就是多项式分解算法，这个可以参考欧几里得多项式除法（https://en.wikipedia.org/wiki/Polynomial_long_division，一个超级叼的算法），总的来说就是算法会不断将原多项式最终分解为下面这样的形式：
(1) 对于最高次为偶数次多项式：
$$ (x^2 + ux + v)(x^2 + u_1x + v_1) ... (x^2 + u_kx + v_k) = 0 $$
(2) 对于最高次为奇数次多项式：
$$ (x^2 + ux + v) ... (x + m) = 0 $$
接下来我们就可以实现算法了，不过我代码整理为基于我自己写得线性代数库了，包含了求逆、向量运算等，然后暂时也没有完全开源，所以，看看算法伪码框架就可以了。
首先是多项式分解算法部分：
b <- edp(u, v, coeff): n <- coeff.size(); b[0 ... n]; rf <- 0; rs <- 1; for i <- 0 ... n: b[i] <- (coeff[n - i - 1] - v * rf) - (u * rs); // 当时是在草稿纸上亲手推的递推式，反正写出来之后会发现有规律，然后就推得这么个算法了 rf <- rs; rs <- b[i]; return b;
然后是寻根部分：
// quads 用于存储分解多项式每个部分的系数部分，最多2个，最少1个 void <- polynomial_product(coeff, max_iter): if coeff.size() < 3: // 如何方程系数个数小于 3 说明是一个一元二次方程，就不用再继续分解了 quads.append(coeff); return; r <- [ 1, 1 ]; epsilon <- 1e-6; i <- 0; while (abs(r[0] + r[1]) > epsilon and i < max_iter): // 多项式分解 a <- edp(u[0], u[1], coeff); b <- [ a.rbegin() + 2 ... a.rend() ]; // 从0到倒数第3个就是分解得到的多项式的系数 // 余项 r <- [ a.end() - 2, a.end() ]; // 剩余的就是余项 cx + d 的系数 c, d if b.size() > 1: z <- [ [-b[0], -b[1]], [0, -b[0]] ]; // 这部分就是 c, d分别对u, v求导得到的矩阵 else: z <- [ [-b[0], -1], [0, -b[0]] ]; // 当分解得到的多项式最高次为 1 时（奇次多项式分解），如 x + m， b_0 = m, b_1 = 1 z <- z.inverse(); // 矩阵求逆 // 牛顿迭代 u <- u - z.dot(r); // 代入迭代公式 i++; quads.append(u); if b.size() > 2: polynomial_product(b); else: quads.append([ b.rbegin() ... b.rend() ]);
最后就只需要遍历 quads 然后对每组使用求根公式得到根即可。