用INTEL CPU计算的结果如下：
Numerator odder N = 9
Denominator Order M = 10
KR =
-118.741401899714 - 141.303692321858i
-118.741401899714 + 141.303692321859i
4655.36084639595 - 1.90177306286648i
4655.36084639598 + 1.90177306348351i
-28916.5722705105 + 18169.1850997184i
-28916.5722705106 - 18169.1850997126i
61276.9997056489 - 95408.5989085931i
61276.9997056478 + 95408.5989085995i
-36902.0468796215 + 196990.463527742i
-36902.0468795457 - 196990.463527748i
ZP =
5.22545336734481 + 15.7295290456402i
5.22545336734481 - 15.7295290456402i
8.77643464007705 + 11.9218538982961i
8.77643464007705 - 11.9218538982961i
10.9343034306161 + 8.40967299601561i
10.9343034306161 - 8.40967299601561i
12.2261314841436 + 5.0127192636517i
12.2261314841436 - 5.0127192636517i
12.8376770778184 + 1.66606258420292i
12.8376770778184 - 1.66606258420292i
显然，后几位有效数字不相同。在几台不同的INTEL CPU上计算得到的结果都相同。但手中只有一台AMD计算机，无法在不同AMD机之间比较，所以也不能结论说AMD浮点运算误差大，而且业界普遍认为AMD浮点运算超强。故敬请有AMD的大虾帮忙难验证一下，下面是MATLAT原码：
function YutouPADE
% =========================================================================
% YUTOUNILT M-file for calulation of Pade approximation
% Note: Pade approximation
%
% exp(z) = Pade(z) = N(z)/M(z) = sum( Ki/(z-pi)
% N(z) = bn*z^n + bn-1*z^n-1 + ... + b1*z + b0
% M(z) = am*z^m + am-1*z^m-1 + ... + a1*z + a0
%
% Last Modified by YuTou v2.5 17-Jul-2007 23:23:28
% =========================================================================
% ----- Direct commands for format dispaly of the results ----
clc
format long g
format compact
% -----
N = input( 'Numerator odder N = ' );
M = input( 'Denominator Order M = ' );
% If you just press the key ENTER, the default value is M=10 and N=9
if isempty(M) || isempty(N)
M = 10;
N = 9;
disp( 'You do not enter any number. System sets M=10 and N=9' );
end
% Coefficients of Polynomial N(z) and M(z)
BN = zeros( 1, N+1 );
AM = zeros( 1, M+1 );
IM = [ M :-1 : 0 ];
IN = [ N :-1 : 0 ];
% Multiplicity 1E-5 is just for accuracy
NF = factorial( N ) .* 1E-5;
MF = factorial( M ) .* 1E-5;
for I = 1 : M+1
TEMP = [ M-IM(I)+1 : M-IM(I)+N ];
AM(I) = PowerNegOne( I ) .* prod( TEMP ) .* 1E-5 .* ( MF./factorial(IM(I)) );
%AM(I) = ((-1)^I) .* factorial( M+N-IM(I) ) .* nchoosek( M, IM(I) );
%AM(I) = ((-1)^I)*factorial(M+N-IM(I))/factorial(M-IM(I))/factorial(IM(I))*MF;
%AM(I) = ((-1)^I)*factorial(M+N-IM(I))*MF/( factorial(M-IM(I))*factorial(IM(I)) );
end
for I = 1 : N+1
TEMP = [ N-IN(I)+1 : N-IN(I)+M ];
BN(I) = prod( TEMP ) .* 1E-5 .* ( NF./factorial(IN(I)) );
%BN(I) = factorial( M+N-IN(I) ) .* nchoosek( N, IN(I) );
%BN(I) = factorial(M+N-IN(I))*NF/( factorial(N-IN(I))*factorial(IN(I)) );
%TempStr = [ TempStr; sprintf('%32.30g', BN(I)) ];
end
% Coefficients of the Partial fractional: KR is resdues and ZP poles
[ KR, ZP, KK ]= residue( BN, AM );
KR
ZP

% --------------------------------------------------
function R = PowerNegOne( I )
if rem( I, 2 ) > 0
R = -1;
else
R = 1;
end