双精度和变精度浮点数测试(Matlab)

                     

贡献者: addis

  • 本词条处于草稿阶段.
预备知识 计算机算数,Matlab 符号计算和变精度计算

代码 1:digits2.m
% 求当前 vpa 变精度计算的有效数字
% digi_dec 是十进制有效数字
% digi_bin 是二进制有效数字
% eps 是最大相对精度
function [digi_dec, digi_bin, eps] = digits2
n = 0;
while true
    n = n + 1;
    if (1 + 2^(-vpa(n))) - 1 == 0
        break;
    end
end
digi_dec = n * log10(2) + 1;
digi_bin = n + 1;
eps = 2^(-n);
end
代码 2:num2bin.m
% x = a*2^Npw (exactly)
% 1 <= a < 2, 精度为 digits2() 输出的 eps
% 若 x 是双精度, 则 a 的绝对精度为 2^(-52), 也就是 eps
% 若 x 是单精度, 则 a 的绝对精度为 2^(-23), 也就是 eps('single')
% 最小的双精度为 2^(-52)*2^(-1022), 最大的双精度为 (2-eps)*2^1023
function [a, Npow] = num2bin(x)
if isnan(x) || isinf(x)
    a = x; Npow = x; return;
end
sgn = sign(x); x = abs(num2sym(x));
Npow = floor(log2(x));
a = sgn * (x * 2^(-Npow));
if a*2^Npow - sgn*x ~= 0
    error('something wrong!');
end
end
代码 3:num2bin2.m
% 精确解出 x = N*2^Npow, 其中 N 是奇数和整数, Npow 是整数
function [N, Npow] = num2bin2(x)
if isnan(x) || isinf(x)
    N = x; Npow = x; return;
end
if x == 0
    N = 0; Npow = 0; return;
end
N = x;
Npow = 0;
for e = [32, 16, 8, 4, 2, 1]
    fac = 2^e;
    while true
        if mod(N, 1) == 0
            break;
        end
        N = N * fac;
        Npow = Npow - e;
    end
    while true
        if mod(N, fac) == 0
            N = N / fac;
            Npow = Npow + e;
        else
            break;
        end
    end
end
end
代码 4:num2sym.m
% 双精度转换为符号表达式
% sym(x) 会对 x 进行猜测和近似,该函数不会
function s = num2sym(x)
[N, Npow] = num2bin2(x);
s = sym(N)*2^sym(Npow);
end
代码 5:num2vpa.m
% 双精度转换为变精度浮点数
% vpa(x) 会对 x 进行猜测和近似,该函数不会
function ret = num2vpa(x, digi)
[N, Npow] = num2bin2(x);
if nargin == 1
    ret = vpa(N)*2^vpa(Npow);
elseif nargin == 2
    ret = vpa(N, digi)*2^vpa(Npow, digi);
end
end

   该函数参考 “计算机算数

代码 6:double2bin.m
% 把一个双精度数转换为 '0' '1' 字符串
% convert a double to bit string
function str = double2bin(x)
if abs(x) < 2^-1022 || isnan(x) || isinf(x)
    error('not implemented!');
end
% 1 bit sign
if x > 0
    str = '0';
else
    str = '1';
end
x = abs(x);
% 11 bit exponent
[~, e] = num2bin(x);
str = [str dec2bin(e + 1023, 11)];
% 52 bit mantessa
i = int64(num2bin2(x));
mstr = dec2bin(i);
str = [str mstr(2:end)];
if numel(str) < 64
    str = [str repmat('0', [1,64-numel(str)])];
elseif numel(str) > 64
    error('something wrong!');
end
end


致读者: 小时百科一直以来坚持所有内容免费,这导致我们处于严重的亏损状态。 长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。 因此,我们请求广大读者热心打赏 ,使网站得以健康发展。 如果看到这条信息的每位读者能慷慨打赏 10 元,我们一个星期内就能脱离亏损, 并保证在接下来的一整年里向所有读者继续免费提供优质内容。 但遗憾的是只有不到 1% 的读者愿意捐款, 他们的付出帮助了 99% 的读者免费获取知识, 我们在此表示感谢。

                     

友情链接: 超理论坛 | ©小时科技 保留一切权利