如何用matlab做高精度计算？【第一辑】

高精度计算是一种程序设计的算法。由于中央处理器的字长限制，如32位CPU中一个整数最大只能取值4,294,967,295（=2^32-1），因此在超范围数值计算中，往往要采用模拟手段。通常使用分离字符的方法来处理数字数组。 维基百科【高精度计算】

对于跟咱一样的普通使用者而言，往往并不关心如何去实现高精度计算，更不会去研究相应的算法。咱这里讲的高精度计算也指的是计算过程中保持数据的精度不丢失。因为内容较多，计划分成三辑进行分享。

第一辑主要介绍matlab自带的高精度计算工具；第二辑主要介绍来自于File Exchange中的两款高精度计算工具箱；第三辑主要介绍一款收费的高精度计算工具箱Multiprecision Computing Toolbox。

matlab自带的高精度计算工具主要依赖于Symbolic Math Toolbox工具箱。因此，想要使用matlab自带的高精度计算工具，务必需要安装Symbolic Math Toolbox工具箱。

两个关键函数 —— vpa与digits。

digits函数使用来控制vpa的计算精度，使用前，按如下方式设置想要的精度即可：digits(num)，num为设置的精度位数。示例如下：

digits(45)  % 
>> vpa(pi)
ans = 3.1415926535897932384626433832795028841971694
>> length('3.1415926535897932384626433832795028841971694')
ans = 45

在计算过程中，精度越高所花费的时间会相应增加，精度越低所花费的时间会相应减少，因此需要在精度与计算时间间做一个权衡。默认情况下，MATLAB使用16位精度。而使用vpa可获得更高的精度，vpa的默认精度为32位，即在未使用digits进行精度位数定义时，digits的返回值为32。如重启matlab后，在命令窗口用π来测试：

>> digits
 Digits = 32
>> dpi = vpa(pi)
dpi = 3.1415926535897932384626433832795
>> length(char(dpi))
ans = 33            % 此处为33是因为有小数点

正如上面所讲matlab的默认浮点数是16位精度，若超过16位则使用四舍五入法进行截断，如下图所示。

那怎么来实现让MATLAB存储高精度数值呢？就得用到咱们主角vpa函数了，使用vpa处理的数据自动为sym型数据。而sym本身也是创建符号变量、表达式、函数、矩阵等函数，且其所创变量、表达式、函数、矩阵同样为sym型数据。

一、创建高精度数据

小伙伴们可能会想直接使用vpa或sym是不是就可以搞定了，那咱们就来看看下面的示例，定义一个高精度小数：3.141514546465512487984541。

>> a = sym(3.141514546465512487984541)
a = 7074061870420537/2251799813685248
>> vpa(a)
 ans = 3.1415145464655123141994863544824
>> b = vpa(3.141514546465512487984541)
b = 3.1415145464655123141994863544824
>> vpa(b)
ans = 3.1415145464655123141994863544824

从结果不难看出，无论直接使用sym还是vpa都未能正确创建咱们想要的数字，这是为什么呢？究其原因，就是MATLAB默认浮点数精度在作怪，超过16位精度的数字就会自动被截断，而其后被无效数值填充，导致无法得出正确结果。那要怎么解决此问题呢？问题根源既然出在数值上，那咱不用数值不就可以了吗？要不试试看？咱直接将数值定义成字符串，然后再使用sym或vpa。

>> a = sym('3.141514546465512487984541')
a = 3.141514546465512487984541
>> b = vpa('3.141514546465512487984541')
b = 3.141514546465512487984541

显然咱的方案是正确的。那有的小伙伴就问了，用num2str来转换上面的数值可以吗？不用怀疑，肯定是不会对的，不信的小伙伴可以私下试试。

二、sym型数据数值化

使用vpa函数可以将sym型数据数值化，通过char函数可将数值化的sym型数据转成可用的字符串型数据。如下：

>> p = sym(pi);
>> a = sym('3.141514546465512487984541');
>> b = sym(1/2);
>> f = a*sin(b*p);
>> val = vpa(f)
val = 3.141514546465512487984541
% vpa有效的保留计算精度，f和a完全相等

三、更改数据的精度

除了使用digits定义精度位数外，使用vpa也可以直接定义计算精度，如取pi的前128位：

>> pi128 = vpa(pi,128)
pi128 = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446

四、符号计算结果数值化

虽然符号计算的结果是精确的，但可能是一种不便于阅读的形式，如使用solve求解高次多项式的根，它给出的不是具体解的值，而是用root来表示根。

>> syms x
y = solve(x^3 - x + 1, x)
y = root(z^3 - z + 1, z, 1)
    root(z^3 - z + 1, z, 2)
    root(z^3 - z + 1, z, 3)
>> vvpa = vpa(y)
vvpa =
0.66235897862237301298045442723905 - 0.56227951206230124389918214490937i
0.66235897862237301298045442723905 + 0.56227951206230124389918214490937i
                                      -1.3247179572447460259609088544781

由于篇幅不宜过长，更多有关vpa的应用小伙伴可以自行研究使用。接下来会在另外两辑中介绍第三方高精度计算工具。欲知后事如何，且看下回分解！

参考资料：

[1] www.mathworks.com/help/symbolic/vpa.html

[2] https://zh.m.wikipedia.org/zh-hans/高精度计算