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

Question

假设您想使用来知道W数字的前有效数字，例如pi vpa。仅仅vpa用那么多的数字打电话是行不通的。考虑以下示例W = 35：

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

四舍五入的原因。具体来说，以上结果似乎表明35pi 的-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

可以看出，

• 从提高所需的小数位数79来80或者81在vpa给出在Matlab相同的答案，因为舍入和运做最后的数字零和Matlab修剪尾随零。
• 八度不修整，因此显示那些零，但它们仍然不正确。

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

---

以上示例说明了

• vpa由于四舍五入，最后一位数字可能不正确；
• 要求多一位数字并不总是足够的。
• 避免舍入问题所需的额外位数可能会很大。当想要的数字后面有很长的九位数时，就会发生这种情况。

那么，是否有一种方法可以W确保数字正确，即不受舍入问题的影响，从而获得数字的前几位有效数字？

Answer 1

首先，它似乎没有可能的情况下，预测vpa将舍远或向零。以下结果与“从零舍入一半”，“从零舍入到偶数”或任何常规规则不一致：

>> 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不想要的部分中初始九的（可能为零）个数，A = 1如果第一个不为9的不想要的数字导致前一个数字舍入为零，A = 0否则为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