matlab - 使用可变精度算法获取数字的第一个有效数字

Question

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

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

这不起作用的原因是舍入。具体来说，上面的结果似乎表明 pi 的第 35 位有效数字是 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

或 Octave :

>> 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

可以看出，

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

因此，在 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