Error with Addition and Subraction of large numbers (e+30)

Question

Ash Ash 2018-3-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/389709-error-with-addition-and-subraction-of-large-numbers-e-30

评论： John D'Errico 2018-3-23

采纳的回答： John D'Errico

在 MATLAB Online 中打开

Hi all Good People of MATLAB

I am facing a problem and I require your help and expertise:

When I perform this computation:

>> 5.357777e30+5.367777e30-5.357777e30-5.367777e30
ans =
    -1.125899906842624e+15

I do not get 0 as I would expect. I was expecting the error for double precision computation to be around 10e-15 to 10e-16 , but apparently representing values in this format (Ae+30) brings up this error.

Would you please help me on what I can do to solve this problem? Is there a way to circumvent this problem and bring the error back to the expected range for double precision without having to resort to vpa() (I am performing the addition and subtraction operations through accumarray() ) or proprietary software like Advanpix?

Or would vpa() be a simple and straightforward solution to this problem? .

Note: The inputs for my computation are not in exponential form "Ae+B" but in double format. However, after several steps of computation, MATLAB automatically converts the values to exponential form when it stores the data into a matrix/variable.

I am facing this problem when I run a matrix of such large numbers through accumarray(), that occurs while I am trying to find the inverse of a 5-dimensional multivector and I am calling the simplifyResult() function from Gerardo Aragon-Camarasa's Clifford Algebra toolbox.

Thank you very much for your help and time.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Guillaume 2018-3-22

在 MATLAB Online 中打开

@Ash Ash, "Note: The inputs for my computation are not in exponential form "Ae+B" but in double format. However, after several steps of computation, MATLAB automatically converts the values to exponential form when it stores the data into a matrix/variable".

matlab does not automatically convert anything, and the way matlab displays the number, in exponential form or otherwise, has no impact whatsoever on how the number is stored in memory and hence never has any impact on the result of any calculation.

Furthermore, regardless of how you enter the number, a double (or single) number is always stored in what you'd probably call an exponential form, in terms of powers of 2. You probably would benefit in learning how double numbers are actually stored in memory.

E.g, the integer 3 is actually stored as

sign: 0
fraction: 2^51
exponent: 1024

which is decoded as

 (-1)^sign * (1 + fraction/2^52) * 2^(exponent-1023) 
 = (-1)^0 * (1 + 0.5) * 2^1
 = 1 * 1.5 * 2

Ash Ash 2018-3-23

I see, thanks for the explanation!

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

John D'Errico 2018-3-21

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/389709-error-with-addition-and-subraction-of-large-numbers-e-30#answer_311119

在 MATLAB Online 中打开

Sorry. You can't do it with a double.

This is a problem known as massive subtractive cancellation. It stems from the fact that a double carries roughly 16 digits in the mantissa. Suppose I try to form the number

x = 1 + 1.001002e-13;
sprintf('%0.55f',x)
ans =
  '1.0000000000001001421168211891199462115764617919921875000'

It really does not matter if the numbers are REALLY large. What matters is what happens in the least significant bits of that number. So if they are all on the order of 1e30, it is the same as if they are all on the order of 1.

So we lost some precision in those TINY bits.

x - 1
ans =
   1.0014e-13

So instead of x-1 being 1.002e-13, the result is not exactly what we want.

The problem arises whenever we subtract two numbers that are very nearly identical, except for their least significant bits.

So, no, you cannot solve the problem using a double.

Worse, you may not even be able to solve it using VPA, UNLESS you are very careful in creating the numbers.

x = sym(1 +  1.0001002e-13)
x =
2251799813685473/2251799813685248
vpa(x-1)
ans =
0.000000000000099920072216264088638126850128174

So vpa got it wrong too! Just throwing vpa at the problem is NOT sufficient.

Suppose I am a bit more careful?

x = sym(1) +  1.0001002e-13
x =
79228162514272261203661565163/79228162514264337593543950336
vpa(x-1)
ans =
0.0000000000001000100200000000039377780293699

So I might have tried this:

x = sym('1 +  1.0001002e-13')
Warning: Support of character vectors will be removed in a future release. Character vectors can be used only for variable names and numbers. Instead, to create symbolic
expressions first create symbolic variables using 'syms'. To evaluate character vectors and strings representing symbolic expressions, use 'str2sym'. 
> In sym>convertExpression (line 1581)
  In sym>convertChar (line 1486)
  In sym>tomupad (line 1236)
  In sym (line 215) 
x =
1.00000000000010001002
vpa(x-1)
ans =
0.00000000000010001001999999999999999999998018

Getting better, except for that irritating warning.

x = 1 + sym('1.0001002e-13')
x =
1.00000000000010001002
vpa(x-1)
ans =
0.00000000000010001001999999999999999999998018

The point is, just throwing vpa at a problem does not always fix it.

I suppose you could use my HPF toolbox too. But even there you need to be careful.

x = hpf(1 + 1.002e-13,50)
x =
1.0000000000001001421168211891199462115764617919922
x-1
ans =
0.00000000000010014211682118911994621157646179199218750000000000

So that got it wrong. But if I create the number in a way that maintains full precision, it works too.

x = hpf('1.002e-13',50) + 1
x =
1.0000000000001002
x - 1
ans =
0.0000000000001002

The point is, you need to understand how numbers are stored in a floating point format, if you are to effectively do computations like this.

11 个评论
显示 9更早的评论隐藏 9更早的评论

John D'Errico 2018-3-21

在 MATLAB Online 中打开

Good question, but it depends on how you create the numbers.

vpi(5.357777e30)
Error using vpi (line 160)
If N is a double, it may be no larger than 2^53 - 1

So, if you create the number as a double, then try to convert it to VPI, it will fail for numbers that large.

Anyway, VPI is pretty old now, and when I wrote it, I did not provide the capability to take a character input for larger numbers. So this fails.

vpi('5.357777e30')
Error using vpi (line 103)
If character supplied, then must be a valid decimal integer

The replacement for VPI is VPIJ (which is fully working as far as I recall, though I did want to rewrite the factor code to use a quadratic sieve.) It can handle string input for large numbers.

vpij('5.357777e30')
ans =
    5357777000000000000000000000000

So, this would work:

vpij('5.357777e30') + vpij('5.367777e30') - vpij('5.357777e30') - vpij('5.367777e30')
ans = 
    0

HPF can work, in several ways.

DefaultNumberOfDigits 100
hpf('5.357777e30') + hpf('5.367777e30') - hpf('5.357777e30') - hpf('5.367777e30')
ans =
    0
sum(hpf([5.357777e30, 5.367777e30, -5.357777e30, - 5.367777e30]))
ans =
    0

John D'Errico 2018-3-21

编辑：John D'Errico 2018-3-21

在 MATLAB Online 中打开

Funny. Looks like I might have a typo in the HPF string parsing code. I'll need to fix that immediately.

Sorry. VPIJ is not on the FEX. I've given it out a few times, but not uploaded it to the FEX. It is quite a bit faster than VPI though.

Regardless, you can just pass a matrix directly to HPF. Converting every number to a string is wrong, and you actually lose accuracy! Num2str is not exact.

A = rand(10);
A_hpf = hpf(A);

hpf here actually converts the numbers exactly as they are stored in the double, as long as you use at least about 55 digits of precision. (I don't recall exactly how many digits are needed.) For example:

format long g
A = rand
A =
       0.744974832237824

But that is NOT the number that is stored in A. We can see that by use of sprintf, with the correct format.

sprintf('%0.55f',A)
ans =
  '0.7449748322378244580477257841266691684722900390625000000'

So now make A an HPF.

DefaultNumberOfDigits 55
hpf(A)
ans =
0.7449748322378244580477257841266691684722900390625000000

Again, num2str is not exact, unless you use the right format spec.

num2str(A)
ans =
    '0.74497'
num2str(A,'%0.55f')
ans =
  '0.7449748322378244580477257841266691684722900390625000000'

Regardless, hpf itself can take the numbers and convert them properly, but as I showed in my answer, you need to be extremely careful in how you convert numbers to ANY high precision form, symbolic or HPF.

Remember that if you create a number as a double, then you are already often lost, because the double has only a limited precision.

I don't have your data matrix, so I cannot comment directly. With HPF, I suggest use of at least 55 digits to be safe, as that always ensures you get the number exactly as it is represented as a double.

John D'Errico 2018-3-22

在 MATLAB Online 中打开

Sometimes an update takes time for others to see. I don't know the magic behind that. I think mirrors might be involved, well, mirror sites.

However, I see an update date of 3/21/18.

It really is not relevant, since you really do not need to convert to strings at all. So that patch will not help you. In fact, that is a bad way to convert a number into HPF form. Just use HPF on the array as you have it as a double. HPF will convert the number into an exact duplicate, without need to use a string.

The only real reason why the string input for HPF should be needed is is you had a long list of digits for some number.

For example, suppose you wanted to form the number:

X = hpf('1.2345678901234567890123456789012')
X =
1.2345678901234567890123456789012

As you can see, the number was entered exactly as we wanted it. However, if you did not use a string form...

X = hpf(1.2345678901234567890123456789012)
X =
1.234567890123456690432135474111419171094894409179687500

The problem is MATLAB sees the number, converts it to a double, then passes it to hpf.

The symbolic toolbox fails in exactly the same way.

X = sym(1.2345678901234567890123456789012)
X =
5559999489923579/4503599627370496
vpa(X)
ans =
1.2345678901234566904321354741114

producing a similarly wrong result.

Ash Ash 2018-3-23

编辑：Ash Ash 2018-3-23

在 MATLAB Online 中打开

I see, thanks for the explanation once again.

John, somehow I am getting this error after downloading the latest version of HPF

hpf('1')
ans =
Error using  - 
Integers can only be combined with integers of the same class, or scalar doubles.
Error in hpf/disp (line 1426)
            F = roundn(F,F.Exponent - F.NumberOfDigits(1));
Error in hpf/display (line 1505)
      disp(F)

========

X = hpf(1.2345678901234567890123456789012)
X =
Error using  - 
Integers can only be combined with integers of the same class, or scalar doubles.
Error in hpf/disp (line 1426)
            F = roundn(F,F.Exponent - F.NumberOfDigits(1));
Error in hpf/display (line 1505)
      disp(F)

John D'Errico 2018-3-23

Huh? How did you install it? None of this makes any sense, unless you did something strange with your search path.

Or, do you have a really old version of MATLAB? It was built in R2012a, at least from what the web page says. I cannot recall when I wrote it. Are you using an older release than that?

请先登录，再进行评论。

Error with Addition and Subraction of large numbers (e+30)

4 个评论
显示 2更早的评论隐藏 2更早的评论

采纳的回答

11 个评论
显示 9更早的评论隐藏 9更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Error with Addition and Subraction of large numbers (e+30)

4 个评论 显示 2更早的评论隐藏 2更早的评论

采纳的回答

11 个评论 显示 9更早的评论隐藏 9更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

4 个评论
显示 2更早的评论隐藏 2更早的评论

11 个评论
显示 9更早的评论隐藏 9更早的评论