谈一谈复数、欧拉公式和割圆术

2020-05-26

引言

写这篇文章的目的是想以一种更简单更易理解也更易实现的方式再阐述一遍经典的割圆术算法,而割圆术算法主要是通过一系列边数越来越多的单位圆内接正n边形来趋近一个圆形,借助于复数和复平面的知识,可以使计算过程中对坐标的处理变得简便.

复数和复平面

在中学时我们学习过了一元二次方程:有一个未知量$x$,已知$a$倍的这个未知量$x$的平方$x^2$加上$b$倍的这个未知量$x$再加上一个常数$c$等于$0$,未知量$x$和$a,b,c,0$这4个数的这种关系

$$ a x^2 + b x + c = 0 $$

提供了关于$x$具体是多少的一条重要的线索,隐式地确定了$x$的.从生活中要考虑的具体的问题可以得到例如

$$ x^2 + 3x + 1 = 0 $$

或者

$$ 10 x + 35 x - 20 = 0 $$

或者

$$ \frac{1}{2} x + 2 x - 1 = 0 $$

这样的具体的方程式,而具体问题的答案则由$x$来确定,而

$$ a x^2 + bx + c = 0 $$

则是这一类方程式的抽象,只要其中的$a,b,c$是实.有时候,对于这样的方程,可以找到实数范围内的一个或者两个$x$使得式子等式左右两端的值相等,例如

$$ x^2 + 3x + 1 = 0 \Longrightarrow x = - \frac{3+\sqrt{5}}{2} \ \text{or} \ x = \frac{\sqrt{5}-3}{2} $$

又例如

$$ 10 x + 35 x - 20 = 0 \Longrightarrow x = -4 \ \text{or} \ x = \frac{1}{2} $$

再例如

$$ \frac{1}{2} x + 2 x - 1 = 0 \Longrightarrow x = -2 - \sqrt{6} \ \text{or} \ x = - 2 + \sqrt{6} $$

我们看到,这些方程式中的$x$都能被解出来,并且这些$x$还都是实数,而每一个实数,都可以在数轴上找到一个点与之对应,如图

数轴

数轴

但有时候,例如当方程为

$$ x^2 + x + 1 = 0 $$

的时候,在实数集里面找不到任何一个实数$x$使等式左右两端的值相等,但是,如果我们想象(image)这样一个数$i$,想象它的平方$i^2$等于$-1$,那么对上述这个方程,我们就能有统一的方法来求解,从而得到

$$ x^2 + x + 1 = 0 \Longrightarrow x = - \frac{1}{2} + \frac{\sqrt{3}}{2} i \ \text{or} \ - \frac{1}{2} - \frac{\sqrt{3}}{2} i $$

这个数$i$是想象出来的,是虚构的,从而被叫做「虚数」(imaginary number),还有无数多个类似$-\frac{1}{2}+\frac{\sqrt{3}}{2}i$或者$-\frac{1}{2}-\frac{\sqrt{3}}{2}i$这样的一元二次方程的解可以通过改变方程$a x^2 + b x + c = 0$中的参数$a,b,c$得到,抽象为

$$ z = p + q i $$

其实就是说,很多$a x^2 + b x + c = 0$的解都具有$z = p + q i$的形式,这里,$i^2 = -1, p, q \in \mathbb{R}$ $z$被叫做「复数」(complex number),因为它是由实数$p$,和虚数$q i$「复合」而成.方程$a x^2 + b x + c = 0$如果有实数解$x$,我们可以把$x$画在一条理论上无限延伸的数轴上,那么,方程$a z^2 + b z + c = 0$如果有复数解$z$,是不是也有一个类似于「数轴」的东西可以把$z$也画出来呢?

是有的.

我们暂时不去讨论为什么复数的数轴不是和实数的数轴一样是一条直线(可以引申出实数是不是和复数一样多的问题),复数的「数轴」叫做「复平面」,类似于中学学到的平面直角坐标系,你知道的,如果一个点$(x_1,y_1)$,那么它可以被画在距离$y$轴的垂直距离为$x_1$,距离$x$轴的垂直距离为$y_1$的地方,如图

笛卡尔坐标系上的一点

笛卡尔坐标系上的一点

非常相似地,如果我们想把一个复数画出来,可以简单地把$z=p+qi$中的$p$,也叫做「实部」(real part),和$q$,也叫「虚部」(imaginary part),分别对应到平面笛卡尔坐标系上的一点$(p,q)$,然后画出来就是

复平面

复平面

这个用来画复数的笛卡尔坐标系不再叫笛卡尔坐标系,而是一般被叫做「复平面」(complex plane),即一个专门用来表述复数的平.这样一来,得到了具体的复数,例如$z_1 = -\frac{1}{2}+\frac{\sqrt{3}}{2}i$,$z_2 = -\frac{1}{2}-\frac{\sqrt{3}}{2}i$,那么如图实数可以在实数轴上被表示出来一样,复数可以在复平面上被表示出来

画在复平面上的复数

画在复平面上的复数

在对复数的几何意义有了一个印象之后,我在此介绍一类特殊的复数,其实就是,形如

$$ z(\theta) = \cos \theta + i \sin \theta $$

这样的复数,这里,自由变化的只有$\theta$,而不再像是在$z = p + i q$中那样让$p$和$q$同时自由变化,所以我们说这是一类特殊的复数,即,诚然,对每一个$\theta \in \mathbb{R}$,我们得到的$\cos \theta + i \sin \theta$都是$p + i q$的形式,这里$p = \cos \theta, q = \sin \theta$,但是,在所有$p \in \mathbb{R}, q \in \mathbb{R}$可能的取值中,不一定都有$p + i q = \cos \theta + i \sin \theta$.

为了加深我们对$\cos \theta + i \sin \theta$这类复数的理解,不妨分别随意取几个$\theta$,例如,分别令$\theta$等于

$$ \frac{1}{12} \pi, \frac{2}{12} \pi, \frac{3}{12} \pi, \cdots, \frac{23}{12} \pi, \frac{24}{12} \pi $$

我们可以得到

$$ \begin{align} & \cos \frac{1}{12} \pi + i \sin \frac{1}{12} \pi \\
& \cos \frac{2}{12} \pi + i \sin \frac{2}{12} \pi \\
& \cos \frac{3}{12} \pi + i \sin \frac{3}{12} \pi \\
& \cdots \\
& \cos \frac{23}{12} \pi + i \sin \frac{23}{12} \pi \\
& \cos \frac{24}{12} \pi + i \sin \frac{24}{12} \pi \end{align} $$

然后再把这些复数,按照之前说的方法,实部对应$x$轴坐标,虚部对应$y$轴坐标,画出来,就可以看到

单位圆上的点

单位圆上的点

可以看到,这些复数画在复平面上,都落在了单位圆上,它们其实都可以看做是复平面上单位圆上的点!

再看,假如我们有一个这样的复数,$z_1 = \cos \theta_1 + i \sin \theta_1$,以及,$z_2 = \cos \theta_2 + i \sin \theta_2$,如果要计算这两个复数的乘积,就可以看到

$$ \begin{align} z_1 z_2 &= (\cos \theta_1 + i \sin \theta_1) (\cos \theta_2 + i \sin \theta_2) \\
&= \cos \theta_1 \cos \theta_2 + \cos \theta_1 i \sin \theta_2 \\
&+ i \sin \theta_1 \cos \theta_2 + i^2 \sin \theta_1 \sin \theta_2 \\
&= \cos \theta_1 \cos \theta_2 + i \cos \theta_1 \sin \theta_2 \\
&+ i \sin \theta_1 \cos \theta_2 + (-1) \sin \theta_1 \sin \theta_2 \\
&= \cos \theta_1 \cos \theta_2 - \sin \theta_1 \sin \theta_2 \\
&+ i(\cos \theta_1 \sin \theta_2 + \sin \theta_1 \cos \theta_2) \\
&= \cos (\theta_1 + \theta_2) + i \sin(\theta_1 + \theta_2) \end{align} $$

$$ (\cos \theta_1 + i \sin \theta_1)(\cos \theta_2 + i \sin \theta_2) = \cos(\theta_1 + \theta_2) + i \sin (\theta_1 + \theta_2) $$

以上这段推导,我们用到的公式有$i^2 = -1$的($i$的定义),还有$\cos (\theta_1 + \theta_2) = \cos \theta_1 \cos \theta_2 - \sin \theta_1 \sin \theta_2$ (和角余弦公式),以及$\sin (\theta_1 + \theta_2) = \sin \theta_1 \cos \theta_2 + \sin \theta_2 \cos \theta_1$(和角正弦公式),这些公式对任意的$\theta \in \mathbb{R}$都是成立的,因而,如果我们记$z(\theta) = \cos \theta + i \sin \theta$,那么$z(\theta_1) z(\theta_2) = z(\theta_1 + \theta_2)$对任意的$\theta_1, \theta_2 \in \mathbb{R}$也都是成立的.

我们可以把$z(\theta) = \cos \theta + i \sin \theta$理解为一根在复平面上指向$\theta$方向的时针,而乘法和除法,则分别是逆时针方向和顺时针方向「拨动」这根时针的一种「运算」,假如说我们有$z_1 = \cos \frac{\pi}{6} + i \sin \frac{\pi}{6}$,我们知道$z_1$指向$\theta = \frac{\pi}{6}$的方向,假如说再有$z_2 = \cos \frac{\pi}{12} + i \sin \frac{\pi}{12}$,并且计算$z_3 = z_1 z_2$,那么我们会得到$z_3 = \cos \frac{\pi}{4} + i \sin \frac{\pi}{4}$,我们知道,按照复平面上复数的画法,$z_3$这根时针是指向$\theta = \frac{\pi}{4}$方向的,等于说,我们把原来指向$30^{\circ}$方向的时针$z_1$,拨动了$z_2$的逆时针$15^{\circ}$,最后得到指向$45^{\circ}$的时针$z_3.推广开来不难理解:既然乘法是逆时针拨动,那么除法肯定就是顺时针拨动了,那么开根号就是半倍角,开三次方根就是让幅角缩小到原来的$\frac{1}{3}$倍,单位长度的复数的乘法和除法无非是幅角的加减运.这里,我们把一个单位长度的复数$\cos \theta + i \sin \theta$的$\theta$叫做幅角,并且这个幅角其实就是复平面上的原点到这个复数$\cos \theta + i \sin \theta$与$1$也就是$\cos 0 + i \sin 0$所夹的角度,所有这类复数$\cos \theta + i \sin \theta$又被叫做单位复数,因为,首先,这些复数到$0$(也就是复平面上的原点)的距离都是$1$,因为一个复数$p + i q$到原点的距离等于

$$ d(p + i q, 0) = \sqrt{p^2 + q^2} $$

所以$\cos \theta + i \sin \theta$到原点的距离就等于

$$ d(\cos \theta + i \sin \theta, 0) = \sqrt{\cos^2 \theta + \sin^2 \theta} = 1 $$

所以被叫做单位复数.

欧拉公式

为了计算单位复数$z(\theta)$的任意$t$次幂$(z(\theta))^t$,我们在此介绍欧拉公式(Euler’s formular),欧拉公式将指数、三角函数和虚数联系起来,对任何一个实数$\theta \in \mathbb{R}$,都有

$$ e^{i \theta} = \cos \theta + i \sin \theta $$

成立.

有了欧拉公式,我们可以将$z(\theta)$展开成

$$ \cos \theta + i \sin \theta $$

再改写成$e^{i \theta}$,其实也就是$z(\theta) = e^{i \theta}.于是乎,$\sqrt{z(\theta)} = \sqrt{e^{i \theta}} = e^{i \frac{\theta}{2}}=z(\frac{\theta}{2})$,也就是

$$ \sqrt{\cos \theta + i \sin \theta} = \cos \frac{\theta}{2} + i \sin \frac{\theta}{2} $$

根据欧拉公式,易知,对任意$t \in \mathbb{R}$,都有$(z(\theta))^t = e^{i t \theta} = z(t \theta)$,直观的理解,一个单位复数的$t$次幂,相当于把它的幅角$\theta$改变到原来的$t$倍,这就是单位复数的幂运算的几何意义.

割圆术

一个单位圆内接正$n$边形是一个多边形,这个多边形的每一个顶点都在同一个单位圆上,并且这个多边形的每一条边都是等长的直线段.

从一个圆内接正$n_0$边形开始($n_0 \geq 3$即可,不管具体是多少),例如说在一次迭代中,记这个多边形的顶点分别为v[1], v[2], ..., v[n0],那么我们可以由v[1], v[2]产生等分圆弧v[1]v[2]的新的顶点v'[1],类似地,再由v[2], v[3]产生v'[2]等分v[2]v[3],一直到产生v'[n0-1]等分v[n0-1]v[n0]v'[n0]等分v[n0]v[1],这样,原来有n0个顶点,又新产生了n0个顶点,我们得到2n0个顶点组成的新的圆内接正多边.

例如,假设最初有v[1], v[2], v[3]圆内接正多边形的顶点把整个圆弧做3等分,那么,经过一次迭代,这三个顶点变为v[1], v'[1], v[2], v'[2], v[3], v'[3],即,总共3个顶点变成总共6个顶点,单位圆内接正3边型变为单位圆内接正6边形,把这个过程重复下去,可以得到单位圆内接正12边型,单位圆内接正24边形,在第$n$次迭代结束后,会得到单位圆内接正$3 \times 2^{n}$边形,并且这个多边形是由顶点表示的.

显然,如下图所示

割圆术逼近一个圆

割圆术逼近一个圆

随着迭代的次数增加,单位圆内接正多边形的边数变得越来越多,这个单位圆内接正多边形看起来也越来越像一个圆,于是,我们可以通过计算最后一次迭代得到的内接正多边形的面积,来估计一个单位圆的面积,而又因为一个单位圆的面积等于$\pi$,从而得以估计$\pi.从上图还可以看到,在第4次迭代时,这个内接正$3 \times 2^4$多边形看起来已经和一个圆没什么分别了.

多亏有了复数的知识和欧拉公式,我们可以很轻松地形式化地描述割圆术涉及到的迭代过程,并且不必去计算复杂的坐标,假设我们从单位圆内接正3边型开始,那么不妨设这三个顶点是

$$ e^{i 0}, e^{i \frac{2 \pi}{3}}, e^{i \frac{4 \pi}{3}} $$

其实刚好就是

$$ 1, \ \cos 120^{\circ} + i \sin 120^{\circ}, \ \cos 240^{\circ} + i \sin 240^{\circ} $$

这三个顶点的绝对位置不重要,重要的仅仅是两两之间到原点的夹角相.我们把当前的顶点列表记做a,列表a的长度记做n,用a[i]表示第i个顶点,接下来我们计算

b[i] = Sqrt[ a[i+1]/a[i] ] * a[i], i = 1,2, ..., n-1
b[n] = Sqrt[ a[n]/a[1] ] * a[n]

这里的Sqrt是开根号函.计算完了之后,可以得到一个新长度同样为n的列表b,对于这个这个新的列表b显然b[i]是等分a[i]a[i+1]表示的圆弧.根据欧拉公式,a[i+1]/a[i]这个新的复数,模长不变(等于1),而幅角等于a[i+1]的幅角减去a[i]的幅角,Sqrt[ a[i+1]/a[i] ]的模长等于a[i+1]/a[i]的模长,而幅角则等于a[i+1]/a[i]的幅角的一半,从而,Sqrt[ a[i+1]/a[i] ]正是圆弧a[i+1]a[i]的等分点.

得到了列表b之后,我们按照

a[1] b[1] a[2] b[2] a[3] b[3] ... a[n-1] b[n-1] a[n] b[n]

的方式合并列表a和列表b,合并后的列表可依同样的规则进行下一次迭代,以此不断重复,直到顶点数足够多(从而边数足够多),而使得误差足够小为止.

例如,进行一次迭代之后

Exp[I * 0], Exp[I * (2 Pi/3)], Exp[I * (4 Pi/3)]

就变为

Exp[I * 0], Exp[I * (Pi/3)], Exp[I * (2 Pi/3)], 
Exp[I * (3 Pi/3)], Exp[I * (4 Pi/3)], Exp[I * (5 Pi/3)]

这里的Exp是指数函数,即$e^x$,而I是单位虚数,即$i, i^2 = -1$.

下面,就此计算方法,我给出具体的计算程序的Mathematica代码:

generate[{prev_, current_}] := {prev, Sqrt[current/prev]*prev, current};

evolve[complexes_] := Flatten[
    Map[
        generate, 
        Transpose[{complexes, RotateLeft[complexes]}]
    ]
] // DeleteDuplicates;

circlePointsInNCuts[nCuts_] := Nest[
    evolve, 
    {
        Exp[I*0], Exp[I*(Pi/2)], Exp[I*Pi], 
        Exp[I*((3/2)*Pi)]
    }, 
    nCuts
];

circlePoints = circlePointsInNCuts[12];

circlePointsToCoords[points_] := (
    {
        Re /@ circlePoints,
        Im /@ circlePoints
    } // Transpose
);

coords = circlePointsToCoords[circlePoints];
circlePolygon = Polygon[coords];

piEstimate = Area[circlePolygon]
errorLog = N[Log[10, Abs[Pi - piEstimate]]]

在实际计算时,从内接正4边形出发,经过12次迭代之后,得到正16384边形,此时的绝对误差约为$10^{-7}$.

误差的估计

因为$\pi$本身不是一个有理数,而我们估计$\pi$实际上是用一个有理数$\hat{\pi}=\frac{p}{q}$去估计,所以,无从得知真实的绝对误差$|\pi - \hat{\pi}| = |\pi - \frac{p}{q}|$具体等于多少.

我们知道,内接多边形当边数足够多时能够近似一个圆,而外接多边形当边数足够多时同样也能够近似一个圆,而外接多边形的面积总是大于圆的面积,圆的面积又总是大于内接多边形的面积,我们记圆的外接n边型为$\Delta_O(n)$,我们记圆为$\bigcirc$,我们记圆的内接n边型为$\Delta_I(n)$,我们记面积函数为$S()$,那么,总是有

$$ S(\Delta_O(n)) \geq S(\bigcirc) \geq S(\Delta_I(n)) $$

并且,当边数$n$趋于无穷时

$$ S(\Delta_O(n)) = S(\bigcirc) = S(\Delta_I(n)) $$

于是,当$n$取值很大很大时,就可以用圆外接多边形的面积来近似圆的面积,即

$$ | S(\bigcirc) - S(\Delta_I(n)) | \approx | S(\Delta_O(n)) - S(\Delta_I(n)) | $$

接下来,让我们来考察一个单位圆的外接正$n$边型,具体地,也是为了简单起见,不妨从$n=3$开始,首先,我们在单位圆上面取3个点,分别是在$v_1 = (1,0)$即$(\cos 0, \sin 0)$,$v_2 = (\cos 120^{\circ}, \sin 120^{\circ})$,和$v_3=(\cos 240^{\circ}, \sin 240^{\circ})$

圆外接正3边形

圆外接正3边形

然后,我们做三条直线,分别经过$v_1, v_2, v_3$,并且分别垂直$o v_1, o v_2, o v_3$, 这三条直线两两之间的交点记做$v1m, v2m, v3m$,易知,$v1m, v2m, v3m$就是外接正3边型的顶点!并且,我们还知道$v1m$和$x$轴的夹角是$\frac{1}{2} \frac{2 \pi}{3}$,但是长度我们待会儿才知道,不妨设为$l.那么$v1m=l(\cos \frac{1}{2} \frac{2 \pi}{3}, \sin \frac{1}{2} \frac{2 \pi}{3})$,我们知道,$v1m - v_2$和$v_2$是垂直的!,那么这两个向量计算向量内积就有

$$ (v1m - v_2) \cdot v_2 = 0 $$

也就是

$$ (l \cos \frac{1}{2} \frac{2 \pi}{3} - \cos \frac{2 \pi}{3}, l \sin \frac{1}{2} \frac{2 \pi}{3} - \sin \frac{2 \pi}{3}) \cdot (\cos \frac{2 \pi}{3}, \sin \frac{2 \pi}{3}) = 0 $$

从这个式子中我们可以解出$l$,最后我们算得$l=2$,也就是说,圆的外接正3边型的顶点到中心的距离是圆的内接正3变形的顶点到中心的距离的2倍!也就是$o v1m = 2 o v_1$.

一般地,对于圆的外接正$n$边型,我们类似地有

$$ (l \cos \frac{1}{2} \frac{2 \pi}{n} - \cos \frac{2 \pi}{n}, l \sin \frac{1}{2} \frac{2 \pi}{n} - \sin \frac{2 \pi}{n}) \cdot (\cos \frac{2 \pi}{n}, \sin \frac{2 \pi}{n}) = 0 $$

可从中解出

$$ l = \frac{(\cos \frac{2 \pi}{n})^2+(\sin \frac{2 \pi}{n})^2}{\cos \frac{\pi}{n} \cos \frac{2 \pi}{n} + \sin \frac{\pi}{n} \sin \frac{2 \pi}{n}} $$

当边数为$n$时,圆外接正$n$边型顶点到圆心的距离为$l$,而当边数为$n$时,圆内接正$n$边型到圆心的距离始终是$1$,我们知道,圆内接正$n$边型和圆外接正$n$边型是相似的,因此,圆外接正$n$边形的面积等于圆内接正$n$边形的面积的$l^2$倍,因此,我们有

$$ S(\Delta_O(n)) = l^2 S(\Delta_I(n)) $$

于是乎

$$ S(\Delta_O(n)) - S(\Delta_I(n)) = (l^2 - 1) S(\Delta_I(n)) $$

于是乎

$$ S(\bigcirc) - S(\Delta_I(n)) \approx S(\Delta_O(n)) - S(\Delta_I(n)) = (l^2 - 1) S(\Delta_I(n)) $$

又因为$S(\bigcirc)=\pi$,所以

$$ \pi - S(\Delta_I(n)) \approx = (l^2 - 1) S(\Delta_I(n)) $$

这里$l$用$n$和三角函数可以计算出来,而$S(\Delta_I(n))$是已经计算出来了的,因而这个误差估计是可计算的.

下图是当边数为$n=4 \times 2^{12} = 16384$时,绝对误差的估计的计算过程

绝对误差的估计

绝对误差的估计

可以看到,用$\pi$的估计值和Mathematica自带的较精确$\pi$直接做比较,绝对误差约为$10^{-7.11348}$,而用我们的算法,绝对误差估计为$10^{-7.43454}$.

完整的代码

下面我将完整的代码呈上

generate[{prev_, current_}] := {prev, Sqrt[current/prev]*prev, current};

evolve[complexes_] := Flatten[
    Map[
        generate, 
        Transpose[{complexes, RotateLeft[complexes]}]
    ]
] // DeleteDuplicates;

circlePointsInNCuts[nCuts_] := Nest[
    evolve, 
    {
        Exp[I*0], Exp[I*(Pi/2)], Exp[I*Pi], 
        Exp[I*((3/2)*Pi)]
    }, nCuts
];

circlePointsToCoords[points_] := {
    Re /@ circlePoints,
    Im /@ circlePoints
} // Transpose

piEstimate[nCuts_ /; nCuts >= 3] := (
    circlePoints = circlePointsInNCuts[nCuts];
    coords = circlePointsToCoords[circlePoints];
    circlePolygon = Polygon[coords];
    N[Area[circlePolygon]]
)

pi = piEstimate[12]

nEdges[nCuts_] := 4*2^nCuts

nEdges[12]

N[Log[10, Abs[Pi - pi]]]

tangentialPolygonRadius[n_ /; n >= 3] := (
    v0 = {1, 0};
    v1 = {Cos[(2 Pi)/n], Sin[(2 Pi)/n]};
    vm = l*{Cos[Pi/n], Sin[Pi/n]};
    lSolution = Solve[(vm - v1).v1 == 0, l][[1]];
    N[l /. lSolution]
)

N[Log[10, (tangentialPolygonRadius[16384]^2) - 1]]

总结

复数和复数的世界奇妙又精彩,我们主要是利用了欧拉公式简化了割圆术迭代过程的表示和计算,最终,得到一系列的顶点,这一系列顶点所描述的多边形和单位圆足够相似,而每一个顶点又可通过复数来描述,也可以轻松地转化成两个实数表示的坐标,此时每个点的两个坐标值都是用三角函数来表示,从而我们将一个圆的面积的估计转换为三角函数函数值的估计的问题,最后,我们巧妙地用外接多边形和内接多边形的差来近似绝对误差,而考虑到外界多边形和内接多边形的面积比刚好等于外接多边形半径的平方,按照这种方法,我们仅需要计算出外接多边形的半径即可,而外界多边形的半径显然是容易计算的,因而误差的估计也变得简单.

数学应用complexeulersformulargeyuanshu

从感知机到神经网络

深入浅出关联规则挖掘和Apriori算法