数学-林士谔算法

代数基本定理

1 代数基本定理

任何复系数一元n次多项式（n至少为1）方程在复数域上至少有一根。

n次复系数多项式方程在复数域内有且只有n个根，重根按重数计算。

证明不会

2 虚根成对定理

在实系数多项式分解中，虚根成对分解，实根单一分解，因此对于奇数次多项式，一定有实根。

简单理解： 假设一个根为虚根 $(x-a-bi)$ ，那么必须有另一个共轭根 $(x-a+bi)$ 这样，两者相乘可以得到实系数二次多项式 $x^2-2ax+a^2+b^2$ ,否则无法消去虚数。

而且值得注意的是任何无实根的实系数多项式 $x^2+px+q$ 都可以分解为$x^2-2ax+a^2+b^2$ . (解出a,b从$q>\frac{p^2}{4}$入手证明)

3 林士谔算法

目标：找到多项式的根（包括实根和虚根）。

由代数基本定理，大于二次的多项式一定存在二次三项式因式。如何在实系数多项式中找到一个根，可以考虑在实系数多项式中找到一个二次三项式的因式。

因此可以用逐次分解二次三项式的方法找到所有根。

原理

$$f(x)=(x^2+px+q)g(x)$$

我们要找到一个二次三项式，直接来找是很难的，可以先试着除一个二次三项式

$$f(x)=(x^2+px+q)g(x)+rx+s$$

我们想办法通过修正p和q来消去余项。

$$p_1=p+dq$$

$$q_1=q+dq$$

r和s由p和q决定，对r和s求偏导，得到关系如下。

$$dr=\frac{\partial r}{\partial p}dp+\frac{\partial r}{\partial q}dq$$

$$ds=\frac{\partial s}{\partial p}dp+\frac{\partial s}{\partial q}dq$$

将f(x)分别对p、q求偏导，得到如下：

$$0=xg(x)+\frac{\partial g(x)}{\partial p}(x^2+px+q)+\frac{\partial r}{\partial p}x+\frac{\partial s}{\partial p}$$

$$0=g(x)+\frac{\partial g(x)}{\partial q}(x^2+px+q)+\frac{\partial r}{\partial q}x+\frac{\partial s}{\partial q}$$

由上二个式子可以变形为

$$xg(x)=-\frac{\partial g(x)}{\partial p}(x^2+px+q)-\frac{\partial r}{\partial p}x-\frac{\partial s}{\partial p}$$

$$g(x)=-\frac{\partial g(x)}{\partial q}(x^2+px+q)-\frac{\partial r}{\partial q}x-\frac{\partial s}{\partial q}$$

可以看到我们要求的偏导数正是xg(x)和g(x)对二次三项式的余式。可以直接计算。

因为我们希望r和s都修正到0，令dr=-r，ds=-s,可以近似方程为

$$-r=\frac{\partial r}{\partial p}dp+\frac{\partial r}{\partial q}dq$$

$$-s=\frac{\partial s}{\partial p}dp+\frac{\partial s}{\partial q}dq$$

可以近似解得dp和dq。

多次近似求解可以不断逼近。

代码如下

#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-6;
double b[2333], c[2333];
void shie(double a[], int n, double &p, double &q)
{
    //b = g(x) = f(x) // (x^2 + px + q)
    memset(b, 0, sizeof(b));
    //c = x^2 * g(x) // (x^2 + px + q) from c[] we can get xg(x) % (x^2 + px + q) and g(x) % (x^2 + px + q)
    memset(c, 0, sizeof(c));
    p = 0, q = 0;
    double dp = 1;
    double dq = 1;
    while (dp > eps || dp < -eps || dq > eps || dq < -eps)
    {
        double p0 = p;
        double q0 = q;
        b[n - 2] = a[n];
        b[n - 3] = a[n - 1] - p0 * b[n - 2];
        c[n - 2] = b[n - 2];
        c[n - 3] = b[n - 3] - p0 * b[n - 2];
        for (int j = n - 4; j >= 0; --j)
        {
            b[j] = a[j + 2] - p0 * b[j + 1] - q0 * b[j + 2];
            c[j] = b[j] - p0 * c[j + 1] - q0 * c[j + 2];
        }
        double r = a[1] - p0 * b[0] - q0 * b[1];
        double s = a[0] - q0 * b[0];
        double rp = c[1];
        double sp = b[0] - q0 * c[2];
        double rq = c[0];
        double sq = -q0 * c[1];
        dp = (rp * s - r * sp) / (rp * sq - rq * sp);
        dq = (r * sq - rq * s) / (rp * sq - rq * sp);
        p += dp;
        q += dq;
    }
}
int main()
{
    double a[] = {-6, 11, -6, 1};
    int n = 3;
    double p, q;
    shie(a, n, p, q);
    cout << p << endl << q;
    return 0;
}

其中计算$x^{2}g(x)$对$(x^2+px+q)$的除法可以得到所有想要结果，非常巧妙。