van Roomen问题
1. 问题
400多年前,数学家Adriaan van Roomen以他的著作Ideae mathematicae中的一个45次方程向当时的欧洲顶尖数学家发起挑战,而挑战的名单中没有法国数学家。荷兰大使据此向法国国王指出:法国没有顶尖数学家,于是法国国王召见François Viète (没错就是韦达定理的那个韦达),让他解这个方程。他靠着窗户站了一会儿后就给出了一个解,并在第二天就给出了另外22个(正)解。
这个方程是这样的[1]:

用现代数学语言表示就是这样一个多项式:
$$P_{45}(x)=45x-3795x^3+95634x^5-1138500x^7+7811375x^9-34512075x^{11}+105306075x^{13}-232676280x^{15}+384942375x^{17}-488494125x^{19}+483841800x^{21}-378658800x^{23}+236030652x^{25}-117679100x^{27}+46955700x^{29}-14945040x^{31}+3764565x^{33}-740259x^{35}+111150x^{37}-12300x^{39}+945x^{41}-45x^{43}+x^{45}$$
方程为:
$$P_{45}(x)=K=\sqrt{\frac{7}{4}-\sqrt{\frac{5}{16}}-\sqrt{\frac{15}{8}-\sqrt{\frac{45}{64}}}}$$
2. 试水
虽然是45次方程,但毕竟是400多年前的。
太简单了!只要运用一下高科技,答案立即就出来了:

Ut legit, ut solvit !! 45个根,23个正根,22个负根,并且有20位的精度。
但是,这毕竟是400多年前的问题,使用计算机数值算法纯属犯规了。虽然没必要像Viète一样用纸笔解方程,还是让我们来推测一下他的解题思路。
首先观察$P_{45}$的这一堆系数。注意到系数的正负号是交替的,那么基本可以断定这个式子是某个函数的展开式,但是究竟是什么函数呢?画个图像看看:

看到这个经典的波形,就知道这一定是某个三角函数了。这样也就能解释为什么这个方程的45个根恰好都是实数了。关键问题就在于如何找到这个三角函数。
先把这一堆系数分解质因数:
$$
\begin{align*}
45&=3^2\times5\\
3795&=3\times5\times11\times23\\
95634&=2\times3^3\times7\times11\times23\\
1138500&=2^2\times3^2\times5^3\times11\times23\\
\cdots&\cdots\\
111150&=2\times3^2\times5^2\times13\times19\\
12300&=2^2\times3\times5^2\times41\\
945&=3^3\times5\times7\\
45&=3^2\times5
\end{align*}
$$
除了质因数大小在一定范围内,似乎没有什么规律。
3. 正式解题
引入三角函数
查阅资料后我发现,这个问题竟与三次方程的三角解法有关。
考察没有二次项的三次方程(Depressed Cubic Equation),即形如$x^3+px+q=0(p,q\not=0)$的方程。Viète对其提出了一种利用三角函数的解法[2]。
关键在于三倍角公式的运用:$4\sin^3\theta=3\sin\theta-\sin3\theta$
令$x=A\sin\theta$,$C=\sin3\theta$,代入三倍角公式:
$$4\left(\frac{x}{A}\right)^3=3\left(\frac{x}{A}\right)-C$$
化简:
$$x^3-\frac{3}{4}A^2x+\frac{1}{4}A^3C=0$$
于是得到:
$$
\begin{cases}
p=-3A^2/4\\
q=A^3C/4
\end{cases}
$$
解这个方程组,得到:
$$
\begin{cases}
A=2\sqrt{-p/3}\\
C=\sin3\theta=4q/A^3=-9q/(2p\sqrt{-3p})
\end{cases}
$$
于是可以求出$\theta$的值[3]:
$$\theta=\frac{1}{3}\arcsin\left(-\frac{9q}{2p\sqrt{-3p}}\right)+\frac{2k\pi}{3},\quad k=0,1,2$$
再代回$x=A\sin\theta$:
$$x=2\sqrt{-\frac{p}{3}}\cdot\sin\left(\frac{1}{3}\arcsin\left(-\frac{9q}{2p\sqrt{-3p}}\right)+\frac{2k\pi}{3}\right),\quad k=0,1,2$$
这个求根公式适用于任意无二次项的三次方程,用公式可求出方程的三个根($k=0,1,2$)。若规定$p<0$,$q<0$,则不等式$-1<-9q/(2p\sqrt{-3p})<1$恒成立,三个根都为实数。
顺便一提,如果将$y=x-a/3$代入一般的三次方程$y^3+ay^2+by+c=0$,即可消去二次项,得到$x^3+(b-a^2/3)x+(2a^3/27-ab/3+c)=0$,因此这个三角函数解法可以直接应用于任意三次方程。
化简多项式
可以使用相同的方法化简$P_{45}$。
既然是45次,那就需要$\sin45\theta$,$45=3^2\times5$,那么显然需要使用一次五倍角公式和两次三倍角公式。
$$\sin45\theta=16\sin^5 9\theta-20\sin^3 9\theta+5\sin 9\theta$$
$$\sin9\theta=-4\sin^3 3\theta+3\sin 3\theta$$
设$a=\sin3\theta$,则:
$$
\begin{aligned}
\sin45\theta&=16(-4a^3+3a)^5-20(-4a^3+3a)^3+5(-4a^3+3a)\\
&=-16384a^{15}+61440a^{13}-92160a^{11}+70400a^9-28800a^7+6048a^5-560a^3+15a
\end{aligned}
$$
设$b=\sin\theta$,则:
$$\sin45\theta=17592186044416 b^{45}-197912092999680 b^{43}+1039038488248320 b^{41}-3380998255411200 b^{39}+7638169839206400 b^{37}-12717552782278656 b^{35}+16168683558666240 b^{33}-16047114509352960 b^{31}+12604574741299200 b^{29}-7897310717542400 b^{27}+3959937231224832 b^{25}-1588210119475200 b^{23}+507344899276800 b^{21}-128055803904000 b^{19}+25227583488000 b^{17}-3812168171520 b^{15}+431333683200 b^{13}-35340364800 b^{11}+1999712000 b^9-72864000 b^7+1530144 b^5-15180 b^3+45 b$$
这个式子竟然比$P_{45}$还要复杂得多,但次数倒确实是45次。
仔细观察便能发现这个多项式的各项系数恰是$P_{45}$系数的$2^n$倍:
$$
\begin{aligned}
45&=45\times2^0\\
15180&=3795\times2^2\\
1530144&=95634\times2^4\\
72864000&=1138500\times2^6\\
\cdots&\cdots\\
3380998255411200&=12300\times2^{38}\\
1039038488248320&=945\times2^{40}\\
197912092999680&=45\times2^{42}\\
17592186044416&=1\times2^{44}
\end{aligned}
$$
设$\sin\theta$展开式的$i$次项系数为$\alpha_i$,$P_{45}$的$i$次项系数为$\beta_i$,则:
$$\alpha_i=\beta_i\cdot2^{i-1}$$
于是容易得到,当$x=2\sin\theta$时:
$$2\sin45\theta=P_{45}(x)$$
于是要解的方程就转化为:$2\sin45\theta=K$
化简常数项
多项式可以化简,这么复杂的常数项一定也可以化简。由于左边是$2\sin45\theta$,可以推测右边也是$2\sin\varphi$形式的。那么
$$\frac{1}{2}K=\frac{1}{2}\sqrt{\frac{7}{4}-\sqrt{\frac{5}{16}}-\sqrt{\frac{15}{8}-\sqrt{\frac{45}{64}}}}=\frac{1}{4}\sqrt{7-\sqrt{5}-\sqrt{30-6\sqrt{5}}}$$
是什么角度的正弦值呢?
注意到(五倍角公式):
$$16\left(\frac{1}{2}K\right)^5-20\left(\frac{1}{2}K\right)^3+5\left(\frac{1}{2}K\right)=\frac{\sqrt{3}}{2}=\sin\frac{\pi}{3}$$
所以:
$$\frac{1}{2}K=\sin\frac{\pi}{15}$$
怎么注意到呢?好吧其实是直接用软件计算$\arcsin K/2$得到12°,再反过来用五倍角公式的。
我在网上找到的所有资料都是计算$\sin12°$,得到其值恰为$K/2$的,似乎没有直接用$K/2$得到12°的方法。说到底这个纯粹是van Roomen为了提高问题难度挖的坑,当时只有对三角函数运用熟练、经验极其丰富的数学家才能够得出结果。
解决问题
两边都化简后,问题就简单了。
$$2\sin45\theta=2\sin\frac{\pi}{15}$$
解得$45\theta=\pi/15+2k\pi$,或$45\theta=14\pi/15+2k\pi$,将将所有满足条件的$\theta$求出,即可得到所有根$x=2\sin\theta$,具体计算就不干了,这种事情还是交给计算机吧。
其实到这里只是使用了多项式关于引入的新变量$\theta$的关系式,而没有得出真正的原函数。通过上面求$x$的过程,很容易得出函数:
$$f(x)=2\sin\left(45\arcsin\left(\frac{x}{2}\right)\right)$$
画出图像:

完美!可以看到函数$f(x)$和$P_{45}$在区间$[-2,2]$是完全重合的。但是$P_{45}$对于$f(x)$不仅是展开式,还是一种解析延拓,在此方程中常数项恰在区间$[-2,2]$之间,否则无法用此方法得到解。
4. 其他方法
复数方法
上面我用倍角公式的方法得到了$P_{45}$的化简,本质上还是误打误撞。事实上,若知道原三角函数,还有一种更直接的使用复数和二项式定理的方法求其展开式,下简述之[2]。
令$t=\arcsin(x/2)$,那么:
$$
\begin{aligned}
\sin\left(45\arcsin\left(\frac{x}{2}\right)\right)&=\sin45t\\
&=\mathrm{Im}(e^{45\mathrm{i}t})=\mathrm{Im}((e^{\mathrm{i}t})^{45})\\
&=\mathrm{Im}((\cos t+\mathrm{i}\sin t)^{45})\\
&=\mathrm{Im}\left(\sum_{k=0}^{45}\binom{45}{k}(\mathrm{i}\sin t)^k\cos^{45-k}t\right)\\
&=\frac{1}{\mathrm{i}}\left(\sum_{k=0}^{22}\binom{45}{2k+1}(\mathrm{i}\sin t)^{2k+1}\cos^{45-2k}t)\right)\\
&=\sum_{k=0}^{22}\binom{45}{2k+1}(-1)^k\sin^{2k+1}t(1-\sin^2t)^{22-k})\\
&=\sum_{k=0}^{22}\binom{45}{2k+1}(-1)^k\sin^{2k+1}t\sum_{j=0}^{22-k}\binom{22-k}{j}(-1)^j\sin^{2j}t\\
&=\sum_{k=0}^{22}\sum_{j=0}^{22-k}\binom{45}{2k+1}\binom{22-k}{j}(-1)^{k+j}\sin^{2k+2j+1}t\\
&=\sum_{k=0}^{22}\sum_{j=k}^{22}\binom{45}{2k+1}\binom{22-k}{j-k}(-1)^{j}\sin^{2j+1}t\\
&=\sum_{j=0}^{22}\sum_{k=0}^{j}\binom{45}{2k+1}\binom{22-k}{j-k}(-1)^{j}\sin^{2j+1}t\\
&=\sum_{j=0}^{22}(-1)^j\left[\sum_{k=0}^{j}\binom{45}{2k+1}\binom{22-k}{j-k}\right]\sin^{2j+1}t
\end{aligned}
$$
所以:
$$
\begin{aligned}
f(x)&=2\sin(45\arcsin(\frac{x}{2}))\\
&=\sum_{j=0}^{22}(\frac{-1}{4})^j\left[\sum_{k=0}^{j}\binom{45}{2k+1}\binom{22-k}{j-k}\right]x^{2j+1}\\
&=P_{45}(x)
\end{aligned}
$$
此方法适用于计算更多类似三角函数的展开式。
递推方法
此外,还有一个相当精彩的方法,通过递推得到了$P_{45}$与三角函数的关系[4]。
定义多项式$P_n$:
$$P_n(x)=\sum_{i=0}^{\lfloor n/2\rfloor}(-1)^i\frac{n}{n-i}\binom{n-i}{i}x^{n-2i}$$
容易知道,原方程左端的多项式恰是这样定义的$P_{45}$。
首先,可以求出$P_n$的递推公式:
$$P_n(x)=xP_{n-1}(x)-P_{n-2}(x)$$
证明如下:
对于最简单的情况,因为:
$$
\begin{aligned}
P_1(x)&=x\\
P_2(x)&=x^2-2\\
P_3(x)&=x^3-3x
\end{aligned}
$$
所以:
$$xP_2(x)-P_1(x)=x(x^2-2)-x=x^3-3x=P_3(x)$$
对于整数$n>3$:
$$
\begin{aligned}
xP_{n-1}(x)&=\sum_{i=0}^{\lfloor(n-1)/2\rfloor}(-1)^i\frac{n-1}{n-1-i}\binom{n-1-i}{i}x^{n-2i}\\
&=x^n+\sum_{i=1}^{\lfloor(n-1)/2\rfloor}(-1)^i\frac{n-1}{n-1-i}\binom{n-1-i}{i}x^{n-2i}\\
P_{n-2}(x)&=\sum_{i=0}^{\lfloor n/2-1\rfloor}(-1)^i\frac{n-2}{n-2-i}\binom{n-2-i}{i}x^{n-2-2i}\\
&=\sum_{i=1}^{\lfloor n/2\rfloor}(-1)^{i-1}\frac{n-2}{n-2-(i-1)}\binom{n-2-(i-1)}{i-1}x^{n-2-2(i-1)}\\
&=\sum_{i=1}^{\lfloor n/2\rfloor}(-1)^{i-1}\frac{n-2}{n-1-i}\binom{n-1-i}{i-1}x^{n-2i}
\end{aligned}
$$
接下来根据奇偶性分类讨论。
若$2\nmid n$,则$\lfloor n/2\rfloor=\lfloor(n-1)/2\rfloor=(n-1)/2$,
$$xP_{n-1}(x)-P_{n-2}(x)=x^n+\sum_{i=1}^{\lfloor n/2\rfloor}(-1)^i\left[\frac{n-1}{n-1-i}\binom{n-1-i}{i}+\frac{n-2}{n-1-i}\binom{n-1-i}{i-1}\right]x^{n-2i}$$
注意中间部分:
$$
\begin{aligned}
&\frac{n-1}{n-1-i}\binom{n-1-i}{i}+\frac{n-2}{n-1-i}\binom{n-1-i}{i-1}\\
=&\frac{n-1}{n-1-i}\cdot\frac{(n-1-i)!}{i!(n-1-2i)!}+\frac{n-2}{n-1-i}\cdot\frac{(n-1-i)!}{(i-1)!(n-2i)!}\\
=&\frac{(n-i)!}{(n-i)i!(n-2i)!}\cdot[\frac{(n-1)(n-2i)}{n-1-i}+\frac{i(n-2)}{n-1-i}]\\
=&\frac{n}{n-i}\cdot\frac{(n-i)!}{i!(n-2i)!}\\
=&\frac{n}{n-i}\binom{n-i}{i}
\end{aligned}
$$
代回原式:
$$
\begin{aligned}
xP_{n-1}(x)-P_{n-2}(x)&=x^n+\sum_{i=1}^{\lfloor n/2\rfloor}(-1)^i\frac{n}{n-i}\binom{n-i}{i}x^{n-2i}\\
&=\sum_{i=0}^{\lfloor n/2\rfloor}(-1)^i\frac{n}{n-i}\binom{n-i}{i}x^{n-2i}\\
&=P_n(x)
\end{aligned}
$$
若$2\mid n$,则$\lfloor n/2\rfloor=\lfloor(n-1)/2\rfloor+1=n/2$,
$$
\begin{aligned}
xP_{n-1}(x)-P_{n-2}(x)=&\sum_{i=0}^{\lfloor n/2\rfloor-1}(-1)^i\frac{n}{n-i}\binom{n-i}{i}x^{n-2i}-(-1)^{n/2-1}\frac{n-2}{n/2-1}\binom{n/2-1}{n/2-1}\\
=&\sum_{i=0}^{\lfloor n/2\rfloor-1}(-1)^i\frac{n}{n-i}\binom{n-i}{i}x^{n-2i}+2(-1)^{n/2}
\end{aligned}
$$
当$i=\lfloor n/2\rfloor=n/2$时,
$$(-1)^i\frac{n}{n-i}\binom{n-i}{i}x^{n-2i}=(-1)^{n/2}\frac{n}{n/2}\binom{n/2}{n/2}x^0=2(-1)^{n/2}$$
所以:
$$xP_{n-1}(x)-P_{n-2}(x)=\sum_{i=0}^{\lfloor n/2\rfloor}(-1)^i\frac{n}{n-i}\binom{n-i}{i}x^{n-2i}=P_n(x)$$
然后来看三角函数中的递推关系。
先从一个最基本的关系开始,即喜闻乐见的积化和差:
$$\cos A\cos B=\frac{\cos(A+B)+\cos(A-B)}{2}$$
整理一下可得:
$$2\cos(A+B)=4\cos A\cos B-2\cos(A-B)$$
令$A=(n-1)\alpha$,$B=\alpha$,则可以得到:
$$2\cos n\alpha=2\cos\alpha\cdot2\cos(n-1)\alpha-2\cos\alpha$$
这时很容易看出此递推关系与$P_n$的相似性,若令$x=2\cos\alpha$,则:
$$
\begin{aligned}
P_1(x)&=x=2\cos\alpha\\
P_2(x)&=x^2-2=4\cos^2\alpha-2=2\cos2\alpha\\
P_3(x)&=x^3-3x=8\cos^3\alpha-6\cos\alpha=2\cos3\alpha
\end{aligned}
$$
对于$n>3$,由于递推关系一致,可以得到如下等式:
$$P_n(2\cos\alpha)=2\cos n\alpha$$
正弦函数的关系式可以用余弦函数的关系式推出。
令$\alpha=\beta-\pi/2$代入上式:
$$
\begin{aligned}
P_n\left(2\cos\left(\beta-\frac{\pi}{2}\right)\right)=&2\cos n\left(\beta-\frac{\pi}{2}\right)\\
P_n(2\sin\beta)=&2\cos\left(n\beta-\frac{n\pi}{2}\right)
\end{aligned}
$$
若$2\nmid n$,则可化为更规整的形式:
$$P_n(2\sin\beta)=2\sin n\beta\cdot(-1)^{\frac{n-1}{2}},\quad2\nmid n$$
事实上,上面的三角函数递推过程就是$n$倍角公式的递推计算过程,这种方法就是将最初的三角函数方法规范化,形成一种确定的技巧。
由这样的推导可知,使用余弦函数也可化简多项式,原函数为$2\cos(45\arccos(x/2))$,而且余弦函数对于此类多项式的使用范围更广,最初使用正弦函数成功化简确实是误打误撞,要不是恰好45满足$2\nmid n$,$2\mid(n-1)/2$,就得不到答案了。(运气真好~)
5. 切比雪夫多项式
在上述递归方法中定义的$P_n$,与第一类切比雪夫多项式很类似,这个多项式通常用$T_n$表示,定义为:
$$
\begin{aligned}
T_0(x)&=1\\
T_1(x)&=x\\
T_n(x)&=2xT_{n-1}(x)-T_{n-2}(x)
\end{aligned}
$$
区别在于递推公式系数中的2,这就导致了两类多项式对应三角函数的区别:
$$
\begin{aligned}
&T_n(x)=\cos(n\arccos(x)),\quad x\in[-1,1]\\
&P_n(x)=2\cos\left(n\arccos\left(\frac{x}{2}\right)\right),\quad x\in[-2,2]\\
&T_n(\cos\theta)=\cos n\theta\\
&P_n(2\cos\theta)=2\cos n\theta
\end{aligned}
$$
切比雪夫多项式于1854年由切比雪夫(Chebyshev)引入,之后在逼近理论中被广泛使用。而究其历史,则早在正式提出的250年前就已被发现与使用,Viète等数学家的洞察力与远见实在令人惊叹。
6. 参考
[1]Adriaan van Roomen. Ideae mathematicae pars prima, sive methodus polygonorum, 1593.
[2]R. L. Herman. van Roomen’s Problem, 2021
http://people.uncw.edu/hermanr/mat346/van_Roomen_s_Problem.pdf
[3]参考文献使用符号$\sin^{-1}$,本文为避免歧义,统一使用符号$\arcsin$
[4]Van Roomen’s Problem: Solution Explained
http://fermatslasttheorem.blogspot.com/2007/02/van-roomens-problem-solution-explained.html