使用可变精度算术获得数字的第一个有效数字

Question

假设您想知道一个数字的第一个W 有效数字，比如pi，使用vpa。简单地用那么多数字调用vpa是行不通的。考虑以下带有W = 35 的示例：

>> disp(vpa(sym('pi'), 35))
3.1415926535897932384626433832795029

这不起作用的原因是四舍五入。具体来说，上面的结果似乎表明 pi 的第 35-th 有效数字是 9，而实际上它是四舍五入的 8：

>> disp(vpa(sym('pi'), 36))
3.14159265358979323846264338327950288

---

从上面看来，一个解决方案似乎是要求一个额外的小数并将其丢弃，这样最后一个幸存的小数就不会出现舍入问题。但这也不在一般情况下起作用，因为四舍五入会导致进位。在 Matlab 中查看这个例子：

>> disp(vpa(sym('pi'), 79))
3.141592653589793238462643383279502884197169399375105820974944592307816406286209
>> disp(vpa(sym('pi'), 80))
3.141592653589793238462643383279502884197169399375105820974944592307816406286209
>> disp(vpa(sym('pi'), 81))
3.141592653589793238462643383279502884197169399375105820974944592307816406286209
>> disp(vpa(sym('pi'), 82))
3.141592653589793238462643383279502884197169399375105820974944592307816406286208999
>> disp(vpa(sym('pi'), 83))
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986

或八度：

>> disp(vpa(sym('pi'), 79))
3.141592653589793238462643383279502884197169399375105820974944592307816406286209
>> disp(vpa(sym('pi'), 80))
3.1415926535897932384626433832795028841971693993751058209749445923078164062862090
>> disp(vpa(sym('pi'), 81))
3.14159265358979323846264338327950288419716939937510582097494459230781640628620900
>> disp(vpa(sym('pi'), 82))
3.141592653589793238462643383279502884197169399375105820974944592307816406286208999
>> disp(vpa(sym('pi'), 83))
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986

可以看出，

• 在 vpa 中将所需小数位数从 79 增加到 80 或 81 在 Matlab 中给出相同的答案，因为四舍五入和进位使最后一位数字为零，而 Matlab 修剪尾随零。
• 八度不修剪，所以它显示那些零，但它们仍然不正确。

因此，无论是在 Matlab 中还是在 Octave 中，正确获取第一个 79 有效数字都需要在这种情况下要求至少 三个额外的数字。

---

以上例子说明了这一点

• vpa的最后几位可能因为四舍五入而丢失；
• 要求一个额外的数字并不总是足够的；
• 避免舍入问题所需的额外位数可以任意大。当想要的数字后面有很长的 9 时，就会发生这种情况。

那么，有没有办法获得一个数字的第一个 W 有效数字，并保证它们是正确的，即不受舍入问题的影响？

Answer 1

首先，似乎无法预测vpa 何时会远离 或趋向 零。以下结果与“离零取整一半”、“取舍一半到偶数”或usual rules 中的任何一个不一致：

>> disp(vpa(sym('0.135'),2))
0.14
>> disp(vpa(sym('0.125'),2))
0.12
>> disp(vpa(sym('0.115'),2))
0.11

Octave 的结果也不一致，和 Matlab 的不同：

>> disp(vpa(sym('0.135'),2))
0.14
>> disp(vpa(sym('0.125'),2))
0.13
>> disp(vpa(sym('0.115'),2))
0.11

这种不可预测性并不会真正影响答案，但它迫使它用更通用的术语来表达，而不假设任何特定的舍入规则。

---

让W 成为想要的位数。超出此范围的所有数字将被称为不需要的。令N 为不需要的部分中初始9 的（可能为零）个数，如果第一个非9 的不需要的数字导致前面的数字从零舍入，则让A = 1，否则为A = 0。该图说明了这一点。

根据问题中的观察，有四种可能的情况。在下面的例子中，想要的位数是W = 3，并且使用了Matlab。

1. N = 0、A = 0：不需要额外的数字。

   >> disp(vpa(sym('0.12345'),3)) % works: first 3 digits are correct
   0.123
   

2. N = 0,A = 1：需要一个额外的数字才能正确获得第一个W = 3数字：

   >> disp(vpa(sym('0.12378'),3)) % doesn't work
   0.124
   >> disp(vpa(sym('0.12378'),4)) % works: first 3 digits are correct
   0.1238
   

3. N > 0, A = 0: N 需要额外的数字：

   >> disp(vpa(sym('0.123994'),3)) % doesn't work
   0.124
   >> disp(vpa(sym('0.123994'),4)) % doesn't work
   0.124
   >> disp(vpa(sym('0.123994'),5)) % works: first 3 digits are correct
   0.12399
   

4. N > 0, A = 1: N+1 需要额外的数字：

   >> disp(vpa(sym('0.123997'),3)) % doesn't work
   0.124
   >> disp(vpa(sym('0.123997'),4)) % doesn't work
   0.124
   >> disp(vpa(sym('0.123997'),5)) % doesn't work
   0.124
   >> disp(vpa(sym('0.123997'),6)) % works: first 3 digits are correct
   0.123997
   

让E 表示需要从vpa 询问的额外数字的数量，以确保第一个W 数字是正确的。那么以上四种情况可以通过E = N + A规则总结。

实际上N 和A 都是未知。因此，一种可能的方法是尝试E = 1 并不断增加E，直到最后获得的数字不是（可能被修剪的）0。然后，丢弃最后的E 数字会得到所需的结果。这种方法使用E = max(1, N+A)；也就是说，除了在实际不需要额外数字时使用一个额外数字（上面的情况 1）之外，额外数字的数量是可能的最小值。

下面的代码为实数或虚数实现了这一点，输出可能使用科学计数法。不支持复数（位数更难从输出字符串中找到）。

s = sym('pi'); % number in symbolic form
W = 79; % number of wanted digits
E = 0; % initiallize
done = false;
while ~done
    E = E+1;
    x = char(vpa(s, W+E));
    y = regexprep(x, '^[+-]?0*|\.0*', ''); % remove sign, leading zeros,
        % decimal point and zeros right after the point; if present
    y = regexprep(y, '\D.*$', ''); % remove exponent and imaginary unit,
        % if present
    num_digits = numel(y); % get number of significant digits in x: 
    done = num_digits==W+E && x(end)~='0'; % the second condition is only
        % required in Octave, but it doesn't harm to keep it in Matlab too
end
c = find(~ismember(x, ['0':'9' '+-.']), 1);
if c % there is an exponent or/and imaginary unit
    result = [x(1:c-1-E) x(c:end)]; % discard last E digits before
        % exponent or imaginary unit
else
    result = x(1:end-E); % discard last E digits
end

Answer 2

出现问题是因为 Matlab 使用guard digits 来提高 Matlab 描述的精度：

您使用 vpa 函数或 digits 函数是保证的位数。在内部， 工具箱可以使用比您指定的更多的数字。这些额外的 数字称为保护数字。

为了解决这个问题，我们必须关闭保护数字的使用。 （不）幸运的是，Matlab 不允许用户指定保护数字的数量，这是内部“计算”的。但是，John D'Errico 写了一个High precision floating point arithmetic class (HPF)，您可以在其中自己指定保护数字的数量。

DefaultNumberOfDigits 100 0
DefaultDecimalBase 1

这会将默认位数设置为 100，保护位数为 0，并将 migits 存储在 1 的“元组”中。如果我们现在回到 Pi 示例，我们得到

pie = hpf('pi',79)
pie = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208

pie = hpf('pi',80)
pie = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089

pie = hpf('pi',81)
pie = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899

pie = hpf('pi',82)
pie = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998 

注意：这是解决舍入问题的解决方案。当您开始使用数字进行计算时，仍然存在标准问题。例如。将数字 5 与 100 数字精度相结合，不会为您提供 sqrt(5) 与 100 数字精度。这基本上就是我们有保护数字的原因，它们在“没有”告诉我们的情况下增加了精度。