matlab 6.5和2010a里的vpa以及pi
最初是matlab中文论坛的网友他寻欢发现的，在matlab 7.0, 7.1, 7.4和2010a下运行vpa(pi,100)分别得到如下结果
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
3.141592653589793115997963468544185161590576171875
2010a的结果跟前几个版本在小数点后15位就不一样了
因为手头没有7.1-7.4的版本，倒是为了用老版的netcdf工具箱还留着个6.5，就在6.5上试了一下，确实跟2010a的结果不一样。很有意思的事情，翻了一下vpa的help，两个版本还是有点区别的。
这是2010a的，6.5的没有红色的那部分
VPA Variable precision arithmetic.
R = VPA(S) numerically evaluates each element of the double matrix
S using variable precision floating point arithmetic with D decimal
digit accuracy, where D is the current setting of DIGITS.
The resulting R is a SYM.
VPA(S,D) uses D digits, instead of the current setting of DIGITS.
D is an integer or the SYM representation of a number.
It is important to avoid the evaluation of an expression using double
precision floating point arithmetic before it is passed to VPA.
For example,
phi = vpa((1+sqrt(5))/2)
first computes a 16-digit approximation to the golden ratio, then
converts that approximation to one with d digits, where d is the current
setting of DIGITS. To get full precision, use unevaluated string or
symbolic arguments,
phi = vpa('(1+sqrt(5))/2')
or
s = sym('sqrt(5)')
phi = vpa((1+s)/2);
Additional examples:
vpa(pi,780) shows six consecutive 9's near digit 770 in the
decimal expansion of pi.
vpa(hilb(2),5) returns
[ 1., .50000]
[.50000, .33333]
See also double, digits.
Overloaded methods:
sym/vpa
从这里来看实际上是vpa(pi,100)这个用法不是很标准，应该是vpa('pi',100)，这样算出来在6.5和2010a下都是
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
到这一步上面的问题基本算结束了，但是为什么2010a下vpa(pi,100)的结果是不正确的，问题出在pi上还是vpa上？
首先想到的是用其他的常数对比一下，比如说自然对数的底e
6.5的
>> vpa(exp(1),100)
ans =
2.718281828459045534884808148490265011787414550781250000000000000000000000000000000000000000000000000
>> vpa('exp(1)',100)
ans =
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427
2010a的
>> vpa(exp(1),100)
ans =
2.71828182845904553488480814849026501178741455078125
>> vpa('exp(1)',100)
ans =
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427
vpa('exp(1)',100)两个是一样的，vpa(exp(1),100)则只差了一串0。后来觉得这个里面还有个常数1，不放心，又试了一下vpa(eps)，两个结果仍然一样，只是显示格式略有区别，而vpa('eps')的结果则是
ans =
eps
被系统按照普通的字符串处理了...
看起来是pi的问题，但pi是内置函数，两个版本pi的help连标点符号都一样。后来不死心，去看vpa的源代码，颇有柳暗花明的感觉
6.5的，真简练
if nargin < 2
d = digits;
end;
r = vpa(sym(s),d);
2010a的
eng = symengine;
if strcmp(eng.kind,'maple')
if nargin == 1
r = mapleengine('vpa',s);
else
r = mapleengine('vpa',s,d);
end
if isa(r,'maplesym')
r = sym(r);
end
else
if nargin < 2
d = digits;
end;
if ischar(s)
ss = evalin(symengine,s);
else
ss = sym(s,'f');
end
r = vpa(ss,d);
end
外层if结构的第一块不用看，红色的部分是关键，当输入是char的时候两个基本的等价(只不过似乎2010a支持多种符号运算引擎，所以多了evalin那行，不过默认安装好像仍然只有一个可以选)，但对于不是字符的情况，6.5下用的是
sym(s)
而2010a下则是
sym(s,'f')
这是sym的help里相关的说明
S = sym(A,flag) converts a numeric scalar or matrix to symbolic form.
The technique for converting floating point numbers is specified by
the optional second argument, which may be 'f', 'r', 'e' or 'd'.
The default is 'r'.
'f' stands for 'floating point'. All values are transformed from
double precision to exact numeric values N*2^e for integers N and e.
'r' stands for 'rational'. Floating point numbers obtained by
evaluating expressions of the form p/q, p*pi/q, sqrt(p), 2^q and 10^q
for modest sized integers p and q are converted to the corresponding
symbolic form. This effectively compensates for the roundoff error
involved in the original evaluation, but may not represent the floating
point value precisely. If no simple rational approximation can be
found, the 'f' form is used.
sym(s)相当于sym(s,'r')这是一种有理数的近似，从help来看，对于可以用两个不太大的p和q以p/q, p*pi/q, sqrt(p), 2^q, 10^q的形式表示的数，r要比f精确，如果找不到这样的形式，则r等价于f，至此，问题全部解决，只是不明白2010a为什么舍弃了精度较高的r而选择f。想要改回来的朋友可以自行修改vpa.m，命令行下open就行了。
个人对符号运算不算熟，有疏漏之处欢迎指正。
最后留个小尾巴，其实6.5和2010a的sym函数也是不一样的
当flag是r的时候没什么区别，f的时候一个用的是2的指数的形式，一个是大整数的分数
6.5
>> sym(1/3,'f')
ans =
'1.5555555555555'*2^(-2)
2010a
>> sym(1/3,'f')
ans =
6004799503160661/18014398509481984
有兴趣可以去研究下~
