数值计算(01):¶

引入:对多项式计算的研究¶

我们先引入我们的研究对象 $$ P(x)=2x^{4}+3x+5x-1 $$}-3x^{2

算法一:暴力计算¶

我们直接就一个一个算那么我们总共要算$4+3+2+1=10$次乘法,$4$次加法这种暴力计算是直白的,但是自然地,我们也能注意到我们在不断地重复做有关$x$幂次的乘法.所以,这一定不是最省力的解法.伴随着多项式阶数$n$的提升,我们将会带来$\frac{n(n-1)}{2}+n$的计算次数,利用渐进分析的知识,我们叫这样的算法是$O(n^2)$的.

算法二:打表法¶

我们注意到,只要计算出$\frac{1}{2}$的一些幂次并打表,我们就能节省一定程度的乘法计算.对于$n$幂次的多项式,我们只需要计算$n-1$次乘法来求得$x$的幂次,并记录.然后求$n-1$个幂次与系数的乘法,最后再算$n-1$次加法就可以了.那么总共就会有$O(3n-3)$的复杂度,根据复杂度分析的理论,我们认为他是一个$O(n)$算法. 这样,我们就将原来的平方复杂度算法转化成了线性复杂度算法.它比较占空间,他的空间和给定的$n$是成正比的.因此出现的下面更加优化的方法.

算法三:秦九韶(Horner)算法¶

以0为基点的秦九韶算法¶

我们必须说明的是,打表法并不是最佳的算法,这是因为打表法需要预留空间来存储计算出的$x$的幂次,这样就会带来更大的空间复杂度. 所谓的Horner方法,就是把已知的多项式写成嵌套乘法的形式. 写法很简单:

取出常数
提一个$x$(未知数)
处理$x$后面的因式
直到提$x$不改变式子结构我们以

$$ P(x)=2x^{4}+3x+5x-1 $$ 为例子: $$ \begin{aligned} 2x}-3x^{2^{4}+3x-3x^{{2}+5x-1&=-1+x(2x}-3x+5)\ &=-1+x(5+x(-3+2x^{2}+3x))\ &=\boxed{-1+x(5+x(-3+x(3+2x)))} \end{aligned} $$ 注意到,这时候我们只要做4次乘法,4次加法就可以实现我们的目标.假如我们的次数提升到}+3x^{2$n$次,那么自然就会有$2n$次运算,那么算法的复杂度就是$O(n)$.但是此时我们的空间复杂度是常数,就是$O(1)$.所以我们在原有的时间复杂度情况下节省了足够的空间.

一般的秦九韶算法¶

在之后的插值计算中,我们有时并不是要大量地提取出$x$,而是可能提取$(x-b_{i})$的这样一个整体.我们叫$b_{i}$为一个基点,$b_{i}[i]$就是我们的基点数组. 当基点数组$b[i]=\{0\}$时,引入基点的Horner形式就会退化成经典Horner形式.

嵌套乘法的MATLAB实现¶

%title:"Extended Horner Method"
%Program 0.1 Nested multiplication
%Evaluates polynomial from nested form using Horner’s Method
%%Input: degree d of polynomial,array of d+1 coefficients c (constant term first),
%x-coordinate x at which to evaluate, and
%array of d base points b, if needed
%Output: value y of polynomial at x
function y=nest(d,c,x,b)
if nargin<4, b=zeros(d,1); end
y=c(d+1);
for i=d:-1:1
y = y.*(x-b(i))+c(i);
end

理想很美好,现实很骨感¶

计算的误差¶

自小便知,差之毫厘,谬之千里.在次数很高的多项式计算中,假如我们输入的是诸如$3$这样的又短又简单的整数,那么我们的算法将完美且高效地输出结果,但是现实中的计算往往更为复杂. 学过计算机基础的同学基本都知道,计算机处理数据存储数据的能力有限,64位计算机接受的无符号整型的最大值是$\text{UINT}\_64=18,446,744,073,709,551,615$,倘若考虑带符号的整型,那么最大只能记录从$-9,223,372,036,854,775,807\sim 9,223,372,036,854,775,807$.如果输入这两个数的绝对值还要大的数,就有可能会溢出,带来不必要的麻烦.同理,计算机能记录的小数位数也是有限的.因此,研究合理表示数值的方法和数值误差的估计就应该提上日程了.

数制¶

我们现在研究合理表示数的方法.大家都知道,现在的数值分析仰赖各种强大的计算软件(以 Matlab,Mathematica,Maple 御三家为首).所以我们的数制要能在计算机上有效实现. 但是,计算机很"机械",说到底,它内部是大量的数字电路,也可以认为就是一堆开关和控制开关的信号.沿用布尔代数的符号逻辑,我们设通为$1$断为$0$,那么计算机能理解的"数"无非是一个01序列.现在,我们的唯一问题是,该怎么样规定这个01序列的规则. 推广地来说,所谓数制,就是为给定的数字符号序列施加一套明确的规则，使之能够唯一地代表某一特定集合或范围内的数,并能进行有效的运算.

十进制:介绍位与权¶

其实介绍到这里,大家应该已经明白我们在数学课上无时无刻不和数制打交道,并且深谙一种数制,那就是十进制(Decimal,DEC).我们学过,个位,十位,百位,千位......因此,这样一个数中单独的数码就叫做位(bit),计算机人一般叫它比特. 而我们刚才讲过,不同的符号通过不同的规则来表示有效的数值,我们学过的个,十,百,千就是这个规则,为权(power).我们的十进制,就是位与位对应权乘积的和,$n$位数写成公式如下

\[ A=\sum_{i=1}^{n} a_{i}10^{i-1} \]

这一点不仅适用于十进整数,又可以扩展用于十进小数,所以任何实数都可以写成这样的形式

\[ A_{Dec}=I+F=\color{red}{d_k 10^k + d_{k-1} 10^{k-1} + \dots + d_1 10^1 + d_0 10^0} +\color{blue}{ d_{-1} 10^{-1} + d_{-2} 10^{-2} + d_{-3} 10^{-3} + \dots} \]

其中$I$一个有限和,而$F$是一个介于$[0,1)$收敛的无穷级数,$d_{i}\in \{0,1,2,3,4,5,6,7,8,9\}$

计算机的思路:二进制¶

现在我们来处理计算机的问题,计算机的数码只有两个,所以首先它的位权得发生改变.由研究十进制的经验我们可以知道: 逢$x$进$1$ $\Leftrightarrow$ 权是$x$的幂次,一般为了区分,我们在二进制前面加$0b$作为区分(计算机人是这样的).由于二进制的权都是$2$的幂次,因而任何二进制实数可以表达成下面的形式 $$ A_{Bin}=I+F=\color{red}{d_k 2^k + d_{k-1} 2^{k-1} + \dots + d_1 2^1 + d_0 2^0} +\color{blue}{ d_{-1} 2^{-1} + d_{-2} 2^{-2} + d_{-3} 2^{-3} + \dots} $$

其中$I$一个有限和,而$F$是一个介于$[0,1)$收敛的无穷级数,$d_{i}\in \{0,1\}$

由上面的研究思路,还可扩展至常用的八进制和十六进制.这不属于数值分析的研究范畴,但是十分有趣.数字电路里我们也常用.

数制转换¶

我们平时常用十进制,电脑常用二进制,如果能有快速转换两者的方法,我们就能更有效地和计算机沟通.

非负十转二:长除法(整数部分)和长乘法(小数部分)¶

假如给定非负数值$135.7$,将它转换成二进制要怎么处理呢

整数部分¶

我们的计算方法就是不断对得到的数求除2的余数(一般称之为取2的模)到为止,然后逆序输出,在这个条件下, 135对应的取法是$1(67),1(33),1(16),0(8),0(4),0(2),0(1),1(0)$,所以我们可以得到: $$ 135 \equiv 0\text{b}10000111 $$ 上面的算法一般称之为长除法.

    graph TD
    A[开始] --> B[/输入非负十进制整数 N/];
    B --> C{N == 0?};
    C -- 是 --> D[输出 0];
    D --> Z[结束];
    C -- 否 --> E[初始化空字符串 BinaryString];
    E --> F{N > 0?};
    F -- 是 --> G[计算 Remainder = N % 2];
    G --> H[BinaryString = Remainder + BinaryString];
    H --> I[N = N / 2];
    I --> F;
    F -- 否 --> J[/输出 "不是非负整数"/];
    J --> Z;

小数部分¶

我们这样计算对应二进制数的小数部分,就是不断地得出的结果乘2,若结果大于1,则记1,否则计0,然后取结果的小数部分继续算

graph TD
    A[开始] --> B[/输入小数部分 F/];
    B --> C[/输入最大迭代次数 MaxIter 或所需精度/];
    C --> D[初始化二进制小数字符串 BinFraction];
    D --> E[初始化计数器 Count 为 0];
    E --> F{F > 0 AND Count < MaxIter?}; 
    F -- 是 --> G[计算 Product = F * 2];
    G --> H{Product >= 1?}; 
    H -- 是 --> I[在 BinFraction 后追加 1];
    I --> J[更新 F = Product - 1]; 
    J --> M[Count = Count + 1];
    H -- 否 --> K[在 BinFraction 后追加 0];
    K --> L[更新 F = Product]; 
    L --> M;
    M --> F;
    F -- 否 --> N[/输出 BinFraction/]; 
    N --> Z[结束];

我们沿着上面的路径来,不考虑MaxIter这样的精度限制,那么就是:$1(1.4),0(0.8),1(1.6),1(1.2),0(0.4),0(0.8),1(1.6),\dots$,所以我们发现它具有循环节,所以 $$ 0.7 \equiv 0\text{b}0.1\overline{0110} $$

结合整数部分的分析,我们有 $$ 135.7 \equiv 0\text{b}10000111.1\overline{0110} $$

非负二转十:直接上定义/解方程¶

有限行二转十其实比十进制转二进制要方便的,其实只要根据二进制数的定义将位权和求出即可.所以我们将重点放在无限循环节的二进小数上,比如$0\text{b}0.1\overline{0110}$,如何消去它的循环节呢?做如下的操作: 我们将循环节乘到前面去,自然地,我们要乘$2^{4}$方.循环节的长度有多长,你就需要乘对应底的$n$次方.进而的得到新的数,以便表达,我们的叫那个循环节$z$. 那么我们就可以列出十进制的方程(就不要列二进制的方程折磨自己了). $$ 16z=z+0\text{b}0110=z+6 $$ 所以可得: $$ z=0.6 $$ 那么怎么和我们原来的数不一样的呢?这是因为我们分析的是没有前置位的循环节,我们实际的做法是 $$ 10110.\overline{0110}-1.\overline{0110}=30z $$ 因此可得: $$ 21=30z $$ 此时便得到正确的答案$0.7$.因此,我们消去循环节的方法就是解方程,假设我们有$m$个前置位$d$和长$n$的循环节$c$,这个部分称为$z$,那么我们就会有 $$ \overline{d_{m}d_{m-1}\dots d_{1}c_{n}c_{n-1}\dots c_{1}}-\overline{d_{m}d_{m-1}\dots d_{1}}=(k^{n+m}-k)z $$ 利用上面的方程,我们就可以解出$z$对应的十进制数.

怎么表达负数? IEEE-754 标准¶

没有规矩,不成方圆.任何形式的有效和高效交流都需要一个固定的标准,因此电气电子工程师学会(IEEE)提出了一种标准:IEEE-754.

科学计数法¶

我们在学数理化的时候,老师经常告诉我们,特别大的数应该用科学计数法来表示,比如我们的 Avogadro 常数 $$ N_{A}=+ 6.02 \times 10^{23} $$ 只要看着这个数值,你就已经知道了 IEEE-754 的三个基础

6.02:一般叫做尾数(Mantissa),一般在1到10之间
23:一般叫做指数(exponent),就是指数函数的那个指数
+:一般叫做符号位(sign),分清数值的正负.

IEEE-754 规定完全符合上面的处理方法.只是由于它服务于电气电子,因而,采用了二进制.

IEEE-754 的规范(选读)¶

现在我们根据上面的认识解释规定

尾数(Mantissa):由于是二进制,因此直接可以写成$1.b_{1}b_{2}...b_{n}...$的形式,其中$b\in \{0,1\}$
指数(exponent):就是一个给定长度的数
符号位(Sign):只有正负,占一个位,0正1负

IEEE-754 规定了不同精度的数的精度情况,下面是 IEEE-754 规定的表格

序号	参数	单精度 (Single)	单精度扩展 (Single Extended)	双精度 (Double)	双精度扩展 (Double Extended)	四倍精度 + (Quadruple +)	扩展精度 # (Extended #)
(1)	$\displaystyle p$ (精度，尾数总位数)	$\displaystyle 24$	$\displaystyle \ge 32$	$\displaystyle 53$	$\displaystyle \ge 64$	$\displaystyle 113$	$\displaystyle 64$
(2)	十进制精度位数 ($\displaystyle p / \log_2(10)$)	$\displaystyle 7.22$	$\displaystyle \ge 9.63$	$\displaystyle 15.95$	$\displaystyle \ge 19.26$	$\displaystyle 34.01$	$\displaystyle 19.26$
(3)	尾数最高有效位 (MSB)	隐藏位 (hidden bit)	未指定 (unspecified)	隐藏位 (hidden bit)	未指定 (unspecified)	隐藏位 (hidden bit)	显式位 (explicit bit)
(4)	实际存储尾数位数	$\displaystyle 23$	$\displaystyle \ge 31$	$\displaystyle 52$	$\displaystyle \ge 63$	$\displaystyle 112$	$\displaystyle 64$
(5)	$\displaystyle E_{max}$ (最大指数)	$\displaystyle +127$	$\displaystyle \ge +1023$	$\displaystyle +1023$	$\displaystyle \ge +16383$	$\displaystyle +16383$	$\displaystyle +16383$
(6)	$\displaystyle E_{min}$ (最小指数)	$\displaystyle -126$	$\displaystyle \le -1022$	$\displaystyle -1022$	$\displaystyle \le -16382$	$\displaystyle -16382$	$\displaystyle -16382$
(7)	指数偏置值 (Exponent bias)	$\displaystyle +127$	未指定 (unspecified)	$\displaystyle +1023$	未指定 (unspecified)	$\displaystyle +16383$	$\displaystyle +16383$
(8)	指数位数	$\displaystyle 8$	$\displaystyle \ge 11$	$\displaystyle 11$	$\displaystyle \ge 15$	$\displaystyle 15$	$\displaystyle 15$
(9)	符号位位数	$\displaystyle 1$	$\displaystyle 1$	$\displaystyle 1$	$\displaystyle 1$	$\displaystyle 1$	$\displaystyle 1$
(10)	格式总位数 ((9) + (8) + (4))	$\displaystyle 32$	$\displaystyle \ge 43$	$\displaystyle 64$	$\displaystyle \ge 79$	$\displaystyle 128$	$\displaystyle 80$
(11)	最大数值范围 (约 $\displaystyle 2^{E_{max}+1}$)	$\displaystyle 3.4028 \times 10^{38}$	$\displaystyle \ge 1.7976 \times 10^{308}$	$\displaystyle 1.7976 \times 10^{308}$	$\displaystyle \ge 1.1897 \times 10^{4932}$	$\displaystyle 1.1897 \times 10^{4932}$	$\displaystyle 1.1897 \times 10^{4932}$
(12)	最小正规格化数值范围 ($\displaystyle 2^{E_{min}}$)	$\displaystyle 1.1754 \times 10^{-38}$	$\displaystyle \le 2.2250 \times 10^{-308}$	$\displaystyle 2.2250 \times 10^{-308}$	$\displaystyle \le 3.3621 \times 10^{-4932}$	$\displaystyle 3.3621 \times 10^{-4932}$	$\displaystyle 3.3621 \times 10^{-4932}$
(13)	最小正非规格化数值范围 (约 $\displaystyle 2^{E_{min}-(\text{行4值})}$)	$\displaystyle 1.4012 \times 10^{-45}$	$\displaystyle \le 1.0361 \times 10^{-317}$	$\displaystyle 4.9406 \times 10^{-324}$	$\displaystyle \le 3.6451 \times 10^{-4951}$	$\displaystyle 6.4751 \times 10^{-4966}$	$\displaystyle 1.8225 \times 10^{-4951}$
(14)	FORTRAN 语言类型	`REAL*4`		`REAL*8`		`REAL*16`	`REAL*10`
(15)	C 语言类型	`float`		`double`		`long double`	`long double`

我们一般而言选用double形.

误差¶

我们的内存不是无穷的,正如你上面的所看到的,我们的 double 型只给我们留了$53$个小数位.对于超过$53$位的小数,我们势必要抛弃一些小数位,进而就造成了误差.

机器精度极限¶

我们在高数阶段学习过$\epsilon-\delta(N)$语言.对于机器来说,它不可能表示无穷小.因此机器判定相等的方式是,只要两个数的差距小于机器精度极限$\epsilon_{mach}$,我们就认为他俩相等了.这是没办法的办法,误差不可避免.机器精度极限是1和离它最近的浮点数的差值的绝对值.对于 double 型而言 $$ \epsilon_{mach}=2^{-52} $$

去尾法¶

有一种很简单的方法是,爷不管了,直接去掉录不下来的位数.这种方法的问题在于,经过去尾的数总是小于/等于原来的数的.倘若对$0.999999999$去尾得到$0$,然而其实际值十分地接近$1$不是嘛?因此,我们需要更加"因材施教"的方法.

(四)舍(五)入法¶

对于十进制数而言,我们见过舍入法,即四舍五入.但是,我们必须强调现在二进制舍入法的一些区别,我们把经过舍入法处理的数$x$记作$fl(x)$

不需要看52位的情况(除了53位还有别的地方为1)¶

此时我们只要看第$53$位的情况,如果第$53$位为$1$,那么对$52$位采用进一法,第$53$位为$0$,那么对$52$位采用去尾法

特殊情况:刚好一半(只有53位为1)¶

在这种出生情况下,我们总是以使得$52$位为偶数的方向舍入,

$52$位为$1$,则向上进一
$52$位为$0$,则向下舍入

方法的优势¶

偏差不会正偏也不会负偏,而是会从统计上相互抵消.

相对误差和绝对误差¶

什么是误差?就是测量值$\text{fl}(x)$和精确值$x_{c}$差的绝对值,这也叫做绝对误差(Absolute Error,AE)

\[ AE=|x_{c}-x| \]

有的时候这两个数都太小,写起来很麻烦,我们一般做个比来体现,称为相对误差(Relative Error) $$ RE=\frac{|\text{fl}(x)-x|}{|x|} $$

近似¶

如果 $$ |\text{fl}(x)-x|\leq \frac{1}{2}\epsilon_{mach} $$

那么叫他近似.

序号	参数	单精度 (Single)	单精度扩展 (Single Extended)	双精度 (Double)	双精度扩展 (Double Extended)	四倍精度 + (Quadruple +)	扩展精度 # (Extended #)
(1)	\(\displaystyle p\) (精度，尾数总位数)	\(\displaystyle 24\)	\(\displaystyle \ge 32\)	\(\displaystyle 53\)	\(\displaystyle \ge 64\)	\(\displaystyle 113\)	\(\displaystyle 64\)
(2)	十进制精度位数 (\(\displaystyle p / \log_2(10)\))	\(\displaystyle 7.22\)	\(\displaystyle \ge 9.63\)	\(\displaystyle 15.95\)	\(\displaystyle \ge 19.26\)	\(\displaystyle 34.01\)	\(\displaystyle 19.26\)
(3)	尾数最高有效位 (MSB)	隐藏位 (hidden bit)	未指定 (unspecified)	隐藏位 (hidden bit)	未指定 (unspecified)	隐藏位 (hidden bit)	显式位 (explicit bit)
(4)	实际存储尾数位数	\(\displaystyle 23\)	\(\displaystyle \ge 31\)	\(\displaystyle 52\)	\(\displaystyle \ge 63\)	\(\displaystyle 112\)	\(\displaystyle 64\)
(5)	\(\displaystyle E_{max}\) (最大指数)	\(\displaystyle +127\)	\(\displaystyle \ge +1023\)	\(\displaystyle +1023\)	\(\displaystyle \ge +16383\)	\(\displaystyle +16383\)	\(\displaystyle +16383\)
(6)	\(\displaystyle E_{min}\) (最小指数)	\(\displaystyle -126\)	\(\displaystyle \le -1022\)	\(\displaystyle -1022\)	\(\displaystyle \le -16382\)	\(\displaystyle -16382\)	\(\displaystyle -16382\)
(7)	指数偏置值 (Exponent bias)	\(\displaystyle +127\)	未指定 (unspecified)	\(\displaystyle +1023\)	未指定 (unspecified)	\(\displaystyle +16383\)	\(\displaystyle +16383\)
(8)	指数位数	\(\displaystyle 8\)	\(\displaystyle \ge 11\)	\(\displaystyle 11\)	\(\displaystyle \ge 15\)	\(\displaystyle 15\)	\(\displaystyle 15\)
(9)	符号位位数	\(\displaystyle 1\)	\(\displaystyle 1\)	\(\displaystyle 1\)	\(\displaystyle 1\)	\(\displaystyle 1\)	\(\displaystyle 1\)
(10)	格式总位数 ((9) + (8) + (4))	\(\displaystyle 32\)	\(\displaystyle \ge 43\)	\(\displaystyle 64\)	\(\displaystyle \ge 79\)	\(\displaystyle 128\)	\(\displaystyle 80\)
(11)	最大数值范围 (约 \(\displaystyle 2^{E_{max}+1}\))	\(\displaystyle 3.4028 \times 10^{38}\)	\(\displaystyle \ge 1.7976 \times 10^{308}\)	\(\displaystyle 1.7976 \times 10^{308}\)	\(\displaystyle \ge 1.1897 \times 10^{4932}\)	\(\displaystyle 1.1897 \times 10^{4932}\)	\(\displaystyle 1.1897 \times 10^{4932}\)
(12)	最小正规格化数值范围 (\(\displaystyle 2^{E_{min}}\))	\(\displaystyle 1.1754 \times 10^{-38}\)	\(\displaystyle \le 2.2250 \times 10^{-308}\)	\(\displaystyle 2.2250 \times 10^{-308}\)	\(\displaystyle \le 3.3621 \times 10^{-4932}\)	\(\displaystyle 3.3621 \times 10^{-4932}\)	\(\displaystyle 3.3621 \times 10^{-4932}\)
(13)	最小正非规格化数值范围 (约 \(\displaystyle 2^{E_{min}-(\text{行4值})}\))	\(\displaystyle 1.4012 \times 10^{-45}\)	\(\displaystyle \le 1.0361 \times 10^{-317}\)	\(\displaystyle 4.9406 \times 10^{-324}\)	\(\displaystyle \le 3.6451 \times 10^{-4951}\)	\(\displaystyle 6.4751 \times 10^{-4966}\)	\(\displaystyle 1.8225 \times 10^{-4951}\)
(14)	FORTRAN 语言类型	`REAL*4`		`REAL*8`		`REAL*16`	`REAL*10`
(15)	C 语言类型	`float`		`double`		`long double`	`long double`