本文记录了两种求解一元三次方程的方法:Cardano (卡丹)公式法和伴随矩阵特征值法。其中详细记录了 Cardano 公式的推导过程,以及使用三角函数对复数进行处理,避免复数运算的方法。然后简要介绍了多项式的伴随矩阵和通过伴随矩阵求特征值来获得根的方法。对各种方法分别给出了使用 C++ 以及 Python 求解单一方程和批量方程的代码。并尝试了使用 Pybind11 将 C++ 代码封装为 Python 模块再调用。最后对不同方法的求解速度进行了对比。文后还有一个彩蛋,介绍了 ChatGPT 给出的求解代码。
通过讨论发现,如果一次性已知大量方程的系数,通过 Cardano 公式使用向量化的 Python 代码求解的速度最快。如果必须通过循环才能获得方程的系数,无法使用向量化操作时,则通过 C++ 实现 Cardano 公式求解单一方程的函数,再封装成 Python 模块后调用的速度最快。本文的所有代码可在文末下载。
Cardano(卡丹)公式
推导过程
任意实系数一元三次方程:
\[ \begin{equation} a x^3 + b x^2 + c x + d = 0, a \neq 0 \end{equation} \]
都可以写为无二次项,三次项系数为 1 的形式:
\[ \begin{equation} y^3 + py + q = 0 \end{equation} \]
其中:
\[ \begin{equation} y = x + \frac{b}{3a}, \ p = \frac{3ac-b^2}{3a^2}, \ q = \frac{2b^3 - 9abc + 27a^2d}{27a^3} \end{equation} \]
Cardano 公式中构造两个新变量 \(u, v\) 使得:
\[ \begin{equation} u + v = y, \ 3uv + p = 0 \end{equation} \]
代入方程可得到
\[ \begin{equation} u^3 + v^3 = -q, \ uv = -\frac{p}{3} \end{equation} \]
可以发现,知道了 \(u^3, v^3\) 的和与积,则他们是下面关于 \(t\) 的二次方程的两个根:
\[ \begin{equation} t^2 + qt - \frac{p^3}{27} = 0 \end{equation} \]
解这个一元二次方程,令
\[ \begin{equation} \Delta = \frac{q^2}{4} + \frac{p^3}{27} \end{equation} \]
即一元二次方程判别式的 4 倍,则根据二次方程的求根公式,有:
\[ \begin{equation} u^3, v^3 = -\frac{q}{2} \pm \sqrt{\Delta} \end{equation} \]
代入原式发现只有当 \( u \neq v \) 时方程成立,则最终得到:
\[ \begin{equation} y = \sqrt[3]{-\frac{q}{2} - \sqrt{\frac{q^2}{4} + \frac{p^3}{27}}} + \sqrt[3]{-\frac{q}{2} + \sqrt{\frac{q^2}{4} + \frac{p^3}{27}}} \end{equation} \]
在 Cardano 公式诞生的时,还没有虚数的概念,因此这一公式就以这种形式表达。
复数域的解
引入复数之后,方程 \(w^3 = 1\) 的解有三个,分别是:
\[ \begin{equation} w_1=1, w_2=\frac{-1+\sqrt{3}i}{2}, w_3=\frac{-1-\sqrt{3}i}{2} \end{equation} \]
它们叫做三次单位根(third root of unity)。我们定义:
\[ \begin{equation} w = \frac{-1+\sqrt{3}i}{2} \end{equation} \]
式 (9) 中加号两边的立方根各有三个根,组合在一起共 9 个根。但是因为 \(u,v\) 的积需要满足式 (5) 中的关系,仅有三个根成立。他们是:
\[ \begin{equation} \begin{aligned} y_1 &= \sqrt[3]{R_1} + \sqrt[3]{R_2} \\ y_2 &= w\sqrt[3]{R_1} + w^2\sqrt[3]{R_2} \\ y_3 &= w^2\sqrt[3]{R_1} + w\sqrt[3]{R_2} \end{aligned} \end{equation} \]
由于 Python 的 numpy
库支持复数运算,所以使用这个公式可以直接在 Python 中获得复数解。
如果不想使用复数运算,则需要对判别式 \(\Delta\) 分类讨论。
当 \(\Delta = 0\) 时,公式可化简为:
\[ \begin{equation} y_1 = 2\sqrt[3]{-\frac{q}{2}}, y_2 = y_3 = -\sqrt[3]{-\frac{q}{2}} \end{equation} \]
方程有三个实根,其中有两个相等。
当 \(\Delta > 0\) 时,\(R_1, R_2\) 都是实数。如果把实部与虚部分开写,并记:
\[ \begin{equation} S = \sqrt[3]{R_1}, T=\sqrt[3]{R_2} \end{equation} \]
则有:
\[ \begin{equation} \begin{aligned} y_1 &= S + T \\ y_2 &= - \frac{1}{2}\left(S + T\right) + \frac{\sqrt{3}}{2}\left(S - T\right)i \\ y_3 &= - \frac{1}{2}\left(S + T\right) - \frac{\sqrt{3}}{2}\left(S - T\right)i \end{aligned} \end{equation} \]
三个解中一个为实数,另外两个为共轭复数。
但如果 \(\Delta < 0\) ,\(R_1, R_2\) 中对负数开平方,无法直接完成实部与虚部的分离。这时,可以借助三角函数来进一步分析。
使用三角函数分离实虚部
当 \( \Delta < 0 \) 时,有 \( p < 0 \),则 \( R_1, R_2\) 是共轭复数。令:
\[ \begin{equation} -\frac{q}{2} \pm i\sqrt{-\Delta} = r(\cos{\varphi} \pm i\sin{\varphi}) \end{equation} \]
根据模相等,再根据实部相等,则有:
\[ \begin{equation} r = \sqrt{-\frac{p^3}{27}}, \ \cos{\varphi} = -\frac{q}{2r} \end{equation} \]
代入 Cardano 公式,得
\[ \begin{equation} y = 2 \sqrt[3]{r} \cos{\frac{\varphi + 2k\pi}{3}}, k=0,1,2 \end{equation} \]
有三个实根。
不仅当 \( \Delta < 0 \) 时可以使用三角函数处理,当 \( \Delta > 0 \) 也可以用以下方法处理。
(1)当 \( \Delta > 0, p < 0 \) 时,有:
\[ \begin{equation} 0 < \sqrt{-\frac{p^3}{27}} < \sqrt{\frac{q^2}{4}} = \frac{|q|}{2} \end{equation} \]
引进辅助角 \(\omega\) :
\[ \begin{equation} \sqrt{-\frac{p^3}{27}} = \frac{q}{2} \sin{\omega}, \ \sqrt{\Delta} = \frac{q}{2}\cos{\omega} \end{equation} \]
则代入
\[ \begin{equation} \begin{aligned} \sqrt[3]{-\frac{q}{2} + \sqrt{\Delta}} = \sqrt[3]{-\frac{q}{2} + \frac{q}{2}\cos{\omega}} = -\sqrt{-\frac{p}{3}} \sqrt[3]{\tan{\frac{\omega}{2}}} \\ \sqrt[3]{-\frac{q}{2} - \sqrt{\Delta}} = \sqrt[3]{-\frac{q}{2} - \frac{q}{2}\cos{\omega}} = -\sqrt{-\frac{p}{3}} \sqrt[3]{\cot{\frac{\omega}{2}}} \end{aligned} \end{equation} \]
再引入角度:
\[ \begin{equation} \tan{\varphi} = \sqrt[3]{\tan{\frac{\omega}{2}}} \end{equation} \]
则有:
\[ \begin{equation} \tan{\varphi} + \cot{\varphi} = \frac{2}{\sin{2\varphi}} \end{equation} \]
令:
\[ \begin{equation} R_0 = -\frac{2\sqrt{-p/3}}{\sin{2\varphi}} \end{equation} \]
得到方程的三个根为:
\[ \begin{equation} \begin{aligned} y_1 &= R_0 \\ y_2 &= -\frac{R_0}{2} + \sqrt{-p}\cot{2\varphi}i \\ y_3 &= -\frac{R_0}{2} - \sqrt{-p}\cot{2\varphi}i \end{aligned} \end{equation} \]
(2)当 \( \Delta > 0, p > 0 \) 时,引进 \(\omega\):
\[ \begin{equation} \sqrt{\frac{p^3}{27}} = \frac{q}{2}\tan{\omega}, \sqrt{\Delta} = \frac{q}{2}\frac{1}{\omega} \end{equation} \]
于是:
\[ \begin{equation} \begin{aligned} \sqrt[3]{-\frac{q}{2} + \sqrt{\Delta}} &= \sqrt[3]{\frac{q\sin^2{\left(\omega / 2\right)}}{\cos{\omega}}} = \sqrt{\frac{p}{3}} \cdot \sqrt[3]{\tan{\frac{\omega}{2}}} \\ \sqrt[3]{-\frac{q}{2} - \sqrt{\Delta}} &= -\sqrt[3]{\frac{q\cos^2{\left(\omega / 2\right)}}{\cos{\omega}}} = -\sqrt{\frac{p}{3}} \cdot \sqrt[3]{\cot{\frac{\omega}{2}}} \end{aligned} \end{equation} \]
再引进角度 \(\varphi\) ,使:
\[ \begin{equation} \tan{\varphi} = \sqrt[3]{\tan{\frac{\omega}{2}}} \end{equation} \]
则有方程的根:
\[ \begin{equation} \begin{aligned} y_1 &= \sqrt{\frac{p}{3}}\left(\tan{\varphi} - \cot{\varphi}\right) = -2\sqrt{\frac{p}{3}}\cot{2\varphi} \\ y_2 &= \sqrt{\frac{p}{3}} \cot{2\varphi} + \frac{\sqrt{p}}{\sin{2\varphi}}i \\ y_3 &= \sqrt{\frac{p}{3}} \cot{2\varphi} - \frac{\sqrt{p}}{\sin{2\varphi}}i \end{aligned} \end{equation} \]
小结
一次三次方程的根的情况为:当 \(\Delta < 0\) 时,有三个实根。当 \(\Delta > 0\) 时,有一个实根和两个共轭虚根。当 \(\Delta = 0\) 时,有三个实根,其中两个相等。
当 \(\Delta > 0\) 时,可以直接完成实虚部分离,也可以借助三角函数完成实虚部分离,前者相对简洁。当 \(\Delta < 0\) 时,需要借助三角函数完成实虚部分离。
综上,把公式总结起来,一元三次方程:
\[ \begin{equation} a x^3 + b x^2 + c x + d = 0, a \neq 0 \end{equation} \]
的解为:
\[ \begin{equation} \begin{aligned} x_0 &= -\frac{b}{3a} \\ p &= \frac{3ac-b^2}{3a^2}, \\ q &= \frac{2b^3 - 9abc + 27a^2d}{27a^3} \\ \Delta &= \frac{q^2}{4} + \frac{p^3}{27} \\ \textrm{If } \Delta \geq 0: \\ S &= \sqrt[3]{-\frac{q}{2} + \sqrt{\Delta}} \\ T &= \sqrt[3]{-\frac{q}{2} - \sqrt{\Delta}} \\ x_1 &= x_0 + S + T \\ x_2 &= x_0 - \frac{1}{2}\left(S + T\right) + \frac{\sqrt{3}}{2}\left(S - T\right)i \\ x_3 &= x_0 - \frac{1}{2}\left(S + T\right) - \frac{\sqrt{3}}{2}\left(S - T\right)i \\ \textrm{Else if } \Delta < 0: \\ r &= \sqrt{-\frac{p^3}{27}} \\ \varphi &= \arccos{\left(-\frac{q}{2r}\right)} \\ x_1 &= x_0 + 2 \sqrt[3]{r} \cos{\frac{\varphi}{3}} \\ x_2 &= x_0 + 2 \sqrt[3]{r} \cos{\frac{\varphi + 2\pi}{3}} \\ x_3 &= x_0 + 2 \sqrt[3]{r} \cos{\frac{\varphi + 4\pi}{3}} \end{aligned} \end{equation} \]
上式有一点需要注意,就是这里的三次根号是在实数范围内的,定义域和值域都是全体实数,且运算结果唯一。在编程时需要对符号进行处理。
代码实现
以下是一个 Cardano 公式的 C++ 函数的实现。将解得的三个根各用实部、虚部两个 double
类型保存,形成一个 6 维数组。
1 | // In file cubicroot_cardano.cpp |
同理,也可以使用以下 Python 代码来实现。源码来自这篇文章 。
1 | import math |
在使用 python 的实现中,直接使用了 complex
数据类型来存储复数,运行函数直接返回三个根。
伴随矩阵特征值
计算原理
对于给定多项式函数:
\[ \begin{equation} P(x) = x^n + c_{n-1}x^{n-1} + c_{n-2}x^{n-2} + … + c_1 x + c_0 \end{equation} \]
其伴随矩阵(Companion Matrix)定义为:
\[ \begin{equation} M_x = \begin{pmatrix} 0 & 0 & \cdots & 0 & -c_0 \\ 1 & 0 & \cdots & 0 & -c_1 \\ 0 & 1 & \cdots & 0 & -c_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \vdots & 1 & -c_{n-1} \end{pmatrix} \end{equation} \]
方程的根就是伴随矩阵的特征值。其简要证明可以在这篇文章找到。这种方式可以求任意次方程的根。
代码实现
求解矩阵的特征值是一个运算量比较大的工作,目前多数库都是调用 Fortran 为底层的 Lapack 工具箱来实现快速计算的。如果使用 C++ 迭代计算速度稍逊,如果使用 Python 迭代计算速度就更不用说了。因此这里不再写具体的矩阵求特征值方法,而是借助 C++ 的 Eigen
模块来完成。在下载之前,需要到其网站下载所需的头文件。在编译 C++ 时,需要使用 -I
标志指示编译器头文件的位置。
1 | // In file: cubicroot_eigen.cpp |
这里使用了与上节相同的函数头,把结果用 6 个值的数组保存。实际上也可以直接返回 Vector3cd
类型的数据。
下面是一段使用 Python 来实现的代码,使用了 numpy
库。
1 | # In file: cubicroot_eigen.py |
批量计算
如果有大量方程需要一次性求解,在 C++ 中可以通过循环完成。而在 Python 中,如果使用 for 循环,效率非常低,需要使用向量化操作,才能大幅提升性能。这里参考这篇文章给出使用 Cardano 公式批量求解的代码:
1 | # In file: cubicroot_cardano_python.py |
这里使用了 mask
来避免用 if
判断,可以提高运算速度。
还有使用特征值法向量化求解的 Python 代码:
1 | # In file: cubicroot_eigen_multi.py |
这里在构建伴随矩阵时,使用了三阶张量。可以一次性求解大批量特征值。
C++ 封装为 Python 模块
使用 Python 与 C++ 对比运行时间不太公平,因为 Python 牺牲了一些速度换来的是极大的便利。下面我们把 C++ 写的函数封装成 Python 模块,再来对比时间。
在 加快 Python 调用 OpenSees 的速度 一文中,我已经介绍了使用原生方法在 Linux 中把 C++ 代码封装成动态链接库,再在 Python 中直接 import
的方法。它通过 Python.h
头文件中定义的 PyModuleDef
PyMethodDef
PyMODINIT_FUNC
等函数和类型来完成封装。
本文介绍另一种方法,使用现在比较流行的 Pybind11 来完成封装工作。
安装 Pybind11
安装的方法有很多,由于我使用的是 Anaconda ,可以直接通过 conda-forge 来安装:
1 | conda install -c conda-forge pybind11 |
此外还需要安装 C++ 的生成工具。如果上文的 C++ 代码可以成功编译,说明已经安装好了 C++ 的生成工具。这里我使用的是 Visual Studio 提供的 MSVC 来构建。这样就不用在 Windows 中再来安装 g++ 了。
MSVC 提供了 cl.exe
link.exe
来完成编译,使用方法与 g++
类似,只是参数的输入方法有所不同。要在终端使用 cl.exe
的话,需要打开 Visual Studio 提供的终端,它会自动设置一些编译所需的环境变量。具体的编译方法可以参考这篇文章。但使用 Pybind11 时不需要我们手动编译,它提供的 extension 可以自动找到系统中的编译器并进行编译,在后文中会提到。但是编写 C++ 函数时需要编译和测试。
写封装用的 C++ 函数
首先创建一个文件夹,命名为 mycubicmodule
,后面把我们编写的模块所需文件都放在这个文件夹里。其结构如下:
1 | mycubicmodule |
其中 cubicroot.cpp
是把前文中 cubicroot_cardan.cpp
这一文件重命名后复制过来的。 bind.cpp
定义了原 C++ 文件与 Python 交互所用的函数,代码如下:
1 | // in mycubicmodule/bind.cpp |
使用 std::vector
类型将返回值封装起来, Pybind11 可以自动识别,并与 Python 中的 list
类型对应。
编写安装文件
在 pyproject.toml
中,根据 PEP 517 的规定,可以定义在使用 pip
安装模块时,安装之前需要加载的模块。这里除了需要使用 setuptools
以外,还需要使用 pybind11
。有了这个文件, pip
在安装模块之前,会自动先安装这两个库,作为我们安装时的依赖项。
1 | [build-system] |
最后是在 setup.py
中定义安装所需的代码:
1 | import sys |
编译安装
以上文件准备好后,即可开始编译。注意计算机中需要有 C++ 编译器,可以使用 Visual Studio 提供的 MSVC 编译器,也可以使用 g++ 等编译器。 Pybind11 会自动在系统中搜索可用的编译器。打开终端,切入 mycubicmodule
文件夹所在的目录,再输入
1 | pip install ./mycubicmodule |
可以看到程序开始自动安装 setuptools
和 pybind11
。安装好后就会使用编译器进行编译,编译的 warning 等信息都会同样打印到控制台。完成之后,输入
1 | pip show mycubicmodule |
即可看到这个库的信息。
运行时间对比
这里对使用 Cardano 公式求解、使用伴随矩阵特征值求解,以及使用 numpy 求解一万个三次方程的速度进行对比。
首先是记录 C++ 求解时间的函数:
1 | // in cubicroot_cardano.cpp and cubicroot_eigen.cpp |
运行 7 次,使用 Cardano 公式求解的平均时间为:1 ms 。使用 Eigen 求解的平均时间为:501 ms 。可见使用 Cardano 公式计算时间非常短。这与 Eigen 求解器使用了比较复杂的数据结构可能也有关系。实际上,使用 clock()
来记录时间的精度不是很高,实际上 Cardano 公式求解的执行时间可能小于 1 ms。
再来看使用 Python 来计算的结果。用 ipython
中的 %%timeit
来记录时间:
1 | # %% |
得到的结果如下:
1 | Time using `numpy.roots`: |
可见最快的是使用向量化 Cardano 公式的 python 函数,平均 1.9 ms 即计算完成。其次是在 C++ 中循环使用 Cardano 公式的函数,平均 5.78 ms 计算完成。再次是 Python 中使用向量化的求特征值方法,平均时间为 10.9 ms 。然后是使用 C++ 单次计算 Cardan 公式,并在 Python 中完成循环,平均时间为 15.6 ms 。其余方法的时间都相对较长。
综上,在一次性知道大批量方程的系数时,使用向量化的 Cardano 公式的运算速度最快。如果一次只知道一个方程的系数,必须通过循环来计算的话,使用 C++ 实现的 Cardano 公式计算比使用 Python 实现的 Cardano 公式计算的效率更高,提升约 7 倍。
彩蛋:咨询 ChatGPT
下面我们来问问 OpenAI 的 ChatGPT ,看他有什么想法,见下图:
可见 ChatGPT 准确地给出了我们刚才讨论的最快方法。把这个函数保存为 chatgpt_suggestion.py
,运行会发现,得出的结果并不正确。也就是说,ChatGPT 依葫芦画瓢给出的代码很像 Cardano 公式,但是实际上并不是准确的 Cardano 公式,不能直接使用。但是它给了一个思路,就是可以直接使用 numpy
的复数运算功能套用 Cardano 公式,而无需分类做各种处理。这是在前文中没有想到的。下面我对 ChatGPT 给出的代码给出评论,再更正为正确的代码。
1 | # In file: chatgpt_suggestion.py |
1 | # %% |
结果为:
1 | Time using ChatGPT suggested direct solver: |
可见向量化运算的速度与我们上文所发现的最快方法所差无几。同时,它直接利用了 numpy
中对负数开根号的运算,代码看起来更简洁。可见 ChatGPT 还是很强大的,只是生成的代码还需要一定的手工更正。
P.S. 也许是运气比较好,第一次问 ChatGPT 就得到了想要的答案。但我发现应该是用 a large number of equations ,而不是用 amount 。但是换用其它词问 ChatGPT 时,很难再得出直接使用 Cardano 公式计算的代码,多数是给使用 np.root
及 scipy.optimize
中的迭代方法等求解。另外,我还尝试了要求 ChatGPT 自行对上面的代码进行更正并加上单元测试,但它似乎没有这个能力。另外 GPT3 还有专门生成代码用的子集 Codex ,目前还不是很好用。我们拭目以待!
以上讨论的代码可以点此下载。
参考文献:
- 华罗庚. 高等数学引论(第一册). 高等教育出版社. 2009.
- Nino Krvavica. On computing roots of quartic and cubic equations in Python. Personal website. 2019. https://nkrvavica.github.io/post/on_computing_roots/
- Mathis Wang. 一元三次方程的求根公式. 知乎专栏:数学学习笔记. 2023. https://zhuanlan.zhihu.com/p/137077558
- Jinyu Li. 多项式的伴随矩阵. Personal website. 2016. https://jinyu.li/2016/12/04/companion-matrix/
- 苏建林. 三次方程的三角函数解法. 科学空间网站. https://www.spaces.ac.cn/archives/831
- Wikipedia. Cubic equation. Wikipedia website. https://en.wikipedia.org/wiki/Cubic_equation
- Weisstein, Eric W. Cubic Formula. From MathWorld–A Wolfram Web Resource. https://mathworld.wolfram.com/CubicFormula.html
- BaseMax. Cubic Equation Calculator. Github repository. https://github.com/BaseMax/CubicEquationCalculator