本文记录了两种求解一元三次方程的方法：Cardano （卡丹）公式法和伴随矩阵特征值法。其中详细记录了 Cardano 公式的推导过程，以及使用三角函数对复数进行处理，避免复数运算的方法。然后简要介绍了多项式的伴随矩阵和通过伴随矩阵求特征值来获得根的方法。对各种方法分别给出了使用 C++ 以及 Python 求解单一方程和批量方程的代码。并尝试了使用 Pybind11 将 C++ 代码封装为 Python 模块再调用。最后对不同方法的求解速度进行了对比。文后还有一个彩蛋，介绍了 ChatGPT 给出的求解代码。

通过讨论发现，如果一次性已知大量方程的系数，通过 Cardano 公式使用向量化的 Python 代码求解的速度最快。如果必须通过循环才能获得方程的系数，无法使用向量化操作时，则通过 C++ 实现 Cardano 公式求解单一方程的函数，再封装成 Python 模块后调用的速度最快。本文的所有代码可在文末下载。

任意实系数一元三次方程：

\[ \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 维数组。

同理，也可以使用以下 Python 代码来实现。源码来自这篇文章 。

这里使用了 mask 来避免用 if 判断，可以提高运算速度。

还有使用特征值法向量化求解的 Python 代码：

这里在构建伴随矩阵时，使用了三阶张量。可以一次性求解大批量特征值。

使用 Python 与 C++ 对比运行时间不太公平，因为 Python 牺牲了一些速度换来的是极大的便利。下面我们把 C++ 写的函数封装成 Python 模块，再来对比时间。

在 加快 Python 调用 OpenSees 的速度 一文中，我已经介绍了使用原生方法在 Linux 中把 C++ 代码封装成动态链接库，再在 Python 中直接 import 的方法。它通过 Python.h 头文件中定义的 PyModuleDef PyMethodDef PyMODINIT_FUNC 等函数和类型来完成封装。

本文介绍另一种方法，使用现在比较流行的 Pybind11 来完成封装工作。

安装的方法有很多，由于我使用的是 Anaconda ，可以直接通过 conda-forge 来安装：

此外还需要安装 C++ 的生成工具。如果上文的 C++ 代码可以成功编译，说明已经安装好了 C++ 的生成工具。这里我使用的是 Visual Studio 提供的 MSVC 来构建。这样就不用在 Windows 中再来安装 g++ 了。

MSVC 提供了 cl.exe link.exe 来完成编译，使用方法与 g++ 类似，只是参数的输入方法有所不同。要在终端使用 cl.exe 的话，需要打开 Visual Studio 提供的终端，它会自动设置一些编译所需的环境变量。具体的编译方法可以参考这篇文章。但使用 Pybind11 时不需要我们手动编译，它提供的 extension 可以自动找到系统中的编译器并进行编译，在后文中会提到。但是编写 C++ 函数时需要编译和测试。

首先创建一个文件夹，命名为 mycubicmodule ，后面把我们编写的模块所需文件都放在这个文件夹里。其结构如下：

其中 cubicroot.cpp 是把前文中 cubicroot_cardan.cpp 这一文件重命名后复制过来的。 bind.cpp 定义了原 C++ 文件与 Python 交互所用的函数，代码如下：

使用 std::vector 类型将返回值封装起来， Pybind11 可以自动识别，并与 Python 中的 list 类型对应。

在 pyproject.toml 中，根据 PEP 517 的规定，可以定义在使用 pip 安装模块时，安装之前需要加载的模块。这里除了需要使用 setuptools 以外，还需要使用 pybind11 。有了这个文件， pip 在安装模块之前，会自动先安装这两个库，作为我们安装时的依赖项。

最后是在 setup.py 中定义安装所需的代码：

以上文件准备好后，即可开始编译。注意计算机中需要有 C++ 编译器，可以使用 Visual Studio 提供的 MSVC 编译器，也可以使用 g++ 等编译器。 Pybind11 会自动在系统中搜索可用的编译器。打开终端，切入 mycubicmodule 文件夹所在的目录，再输入

可以看到程序开始自动安装 setuptools 和 pybind11 。安装好后就会使用编译器进行编译，编译的 warning 等信息都会同样打印到控制台。完成之后，输入

即可看到这个库的信息。

这里对使用 Cardano 公式求解、使用伴随矩阵特征值求解，以及使用 numpy 求解一万个三次方程的速度进行对比。

首先是记录 C++ 求解时间的函数：