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

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.

[Stephen23]

It seems like quite a reasonable magnitude calculation error to me:

>> eps(5.357777e30)
ans =   1.1259e+015

"I do not get 0 as I would expect"

Why do you expect to get zero when you are using double values well above the flintmax limit?

[Ash Ash]

I'm sorry, I did not know about flintmax. Thank you for pointing this out to me. When I wrote "expect" I was referring to the expected mathematical result on pen and paper.

[Guillaume]

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

I see, thanks for the explanation!

Answer 1

[John D'Errico]

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.

[Stephen23]

The values seem to be integer: would vpi be an option?

[John D'Errico]

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

[Ash Ash]

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

Wow, Thank you for your long detailed explanation. Honestly, I am new to this problem.

What is your recommendation on the best way to work with this problem in accumarray()? Can I convert the input array (eg 100x1 matrix) into a vpij or hpf data format before passing it into accumarray? Do you recommend working with vpij or HPF?

PS: I did a search for 'vpij' but I was not able to find any results at FileExchange.

Thank you.

[Ash Ash]

Hi, I was trying to work with HPF, by converting a double matrix to a hpf matrix as an input to accumarray:

for i=1:size_in
   in_hpf(i)=hpf(num2str(in(i,1)));
end
=====================================================
Error using find
Second argument must be a positive scalar integer.

Error in hpf (line 2997)
                k = find(D ~= 0,'1','last');

Would you please advise me on how I can properly do this?

This code without the num2str conversion does not produce an error but it is not able to produce my desired result with accumarray.

for i=1:size_in
    in_hpf(i)=hpf(in(i,1));
end

Edit: I can't even convert a single character array into hpf. Am I missing something? Would you please guide me?

>> hpf('5e22')
Error using find
Second argument must be a positive scalar integer.

Error in hpf (line 2997)
                k = find(D ~= 0,'1','last');

Thank you.

[Steven Lord]

The second input to find is supposed to be a number, not the text representation of a number.

x = 1:10;
x(x < 5) = 0;
find(x, 1, 'first')

[John D'Errico]

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]

New version of HPF posted.

[Ash Ash]

Thanks John, I went to the HPF FEX page here but it says that it was updated on 11 Nov 2015, and the problematic line is still there:

   k = find(D ~= 0,'1','last');

[John D'Errico]

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]

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]

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?