How to sum up large values with high precision?

Hello all together,

I would like to calculate a formula which consists of the alternating (positive/negative) summation of very large values. For small values of N (for example N=3) ist works well. Also for larger values of N, if the gaps between the "Mueh" values are large enough. For the planned application of the formula, however, the values must be quite close together and N should nevertheless be able to reach at least values of 20 or 30. Here´s my current version using symbolic functions and vpa:

% Definition of values
N=11;
Mueh=[0.83 0.84 0.85 0.86 0.87 0.89 0.90 0.91 0.92 0.93 0.94];
Muehc=num2cell(Mueh);
T=25;

% Creating symbolic variables
syms C zaehler nenner CfF MfF;
mueh=sym('m',[1,N]);
syms t;

% Symbolic computing of the C-Values
for k=1:1:N
    C(1,k)=1;
end
for i=2:1:N
    for k=1:1:N
        zaehler=1;
        nenner=1;
        for j=0:1:(i-2)
            zaehler=vpa((zaehler*mueh(N-j)),50);
        end
        for l=0:1:(i-1)
            if (l~=(k-1))
                nenner=vpa((nenner*(mueh(N-l)-mueh(N-k+1))),50);
            end
        end
        C(i,k)=vpa((zaehler/nenner),50);
    end
end

% Symbolic computing of the function (here only for n=1)
for n=1:1:1

      for j=1:1:N-n+1
          CfF(n,j)=vpa(C((N+1-n),j),50);
      end

      for l=1:1:N-n+1
          MfF(n,l)= vpa(mueh(N-l+1),50);
      end

    f(t,mueh)=sum(vpa((vpa(CfF(n,:),50)) .* vpa(exp(-vpa(MfF(n,:),50)*t),50),50))
    h(t)=f(t,Muehc{:}) 
  end

% Plot h

fplot(h, [0 5]);

The current result - although using sym and vpa looks like this:

I also tryed XSum without success.

Do you have an idea how to fix it?

Many thanks and best regards, Johannes

Answer 1

[John D'Errico]

There are a LOT of problems in this code, if you are trying to obtain high precision in the result.

To start:

Mueh=[0.83 0.84 0.85 0.86 0.87 0.89 0.90 0.91 0.92 0.93 0.94];

You may THINK you have created those values exactly. But you did not. In fact, they are approximate, correct to roughly 16 decimal digits.

sprintf('%0.55f',Mueh(1))
ans =
    '0.8299999999999999600319711134943645447492599487304687500'

As you can see, what MATLAB stored was only an approximation. The same is true for the rest of that vector. You can do further computations with those numbers, but if you start with a double precision approximation, you will never have more than the 16 significant digits you stated with. A 50 digit computation is just then a waste of time and CPU cycles.

Next, I see you doing things like this:

% Creating symbolic variables
syms C zaehler nenner CfF MfF;

But then you assign them with values as

zaehler=1;

WRONG. You are thinking as if that variable is forever after defined as symbolic! Not true. Using syms is NOT a permanent declaration of class, something that sticks forever for that variable.

clear
syms X
whos X
Name      Size            Bytes  Class    Attributes
X         1x1                 8  sym

So MATLAB knows that X is symbolic. However, if I now assign a double into X,

X = 3;

whos X
Name      Size            Bytes  Class     Attributes
X         1x1                 8  double

X is now a double precision number.

When you make an assignment like that, X gets overwritten. Everything about X gets overwritten.

I also see you trying to index into an array using a symbolic index, t. You can't do that.

Let me next see if I can re-write your code so it will be more accurate...

I'll use my own hpf tools.

% Definition of values
DefaultNumberOfDigits 100 5
N=11;
mueh=hpf(83:94,100)/100; % These are now exact.

% Compute C in HPF form
C = zeros(N,N,'hpf');
C(1,:) = 1; % small integers are exact, as flints

for i=2:1:N
  for k=1:1:N
      zaehler=hpf(1);
      nenner=hpf(1);
      for j=0:1:(i-2)
          zaehler=zaehler*mueh(N-j);
      end
      for l=0:1:(i-1)
          if (l~=(k-1))
              nenner=nenner*(mueh(N-l)-mueh(N-k+1));
          end
      end
      C(i,k)=zaehler/nenner;
  end
end

So the above computation of C is done correctly, carrying 100 digits, with 5 guard digits as insurance. That is what the statement

DefaultNumberOfDigits 100 5

does for you. The 5 is the number of extra guard digits, so you could specify more if that would be useful.

You should be able to do the remainder.

[Johannes Uhl]

Thank you very much for your helpful answer! I followed your instructions and wrote the last formula using hpf and it works perfect!

Answer 2

[Michael Doxey]

Hi Johannes,

I understand that you are trying to create a formula adding large values, but you are currently running into an issue where the results are not accurate for these large values. I have a couple suggestions that may be able to help you with getting higher precision when summing large values. Here are a few different approaches that you may try:

1. Use the 'digits' function, which can be found here. For example:

>> digits(100);
>> a=vpa('115792089237316195423570985008687907853269984665640564039457584007913129639936');
>> b=vpa(a+1,100);
>> b-a
ans = 
1.0

2. Use "symbolic numbers", which can be specified using character vectors, such as:

>> a = sym('10000000000000000000000000000000000000000000000000');
>> b = a+1;
>> b-a
ans =
1

3. Make use of the File Exchange submission Variable Precision Integer Arithmetic, found here.

EDIT: As John D'Errico kindly pointed out below, it seems you are working with floating point numbers, so it seems that the HPF class may be of more use to you.

Hope these suggestions can help increase your precision when summing large values.

[Karan Gill]

Great answer! Just a quick note that "digits", "vpa", and "sym" are all symbolic functions using symbolic values. However, "vpa" is meant for high-precision floating-point arithmetic, while "sym" is meant for exact symbolic values. Choose one of them based on what you want.

[John D'Errico]

VPI is for large integer arithmetic only, but it appears this is a floating point question. For high (user specified) precision floating point numbers, you can use HPF.

[Walter Roberson]

'However, "vpa" is meant for high-precision floating-point arithmetic, while "sym" is meant for exact symbolic values'

No, sym() is the underlying objects used by vpa() . sym() has provision for converting MATLAB numeric values of all numeric classes into symbolic numbers, with there being options to control how the conversion is done. Also sym() can take a character vector representing a number and will convert quoted numbers into symbolic floating point if appropriate. sym('1.37') for example.

[Karan Gill]

@Walter, I'm afraid I didn't understand how your explanation conflicts with mine. I agree with what you said. As for my reply, I was pointing this out. Or did I get something wrong?

>> sym(1.37) % exact
ans =
137/100
>> vpa(1.37) % floating-point
ans =
1.37

[Walter Roberson]

It is difficult to say that sym(1.37) is exact. sym() with the default 'r' conversion looks for a fraction that is "sufficiently close" to the input value, with the factions being examined including multiples of Pi and including square roots.

The output value from the conversion is an extended rational (extended because pi and square roots are not rationals) that exactly represent the approximation -- but it is an approximation. For example,

>> sym(3^(1/3))
ans =
3247657313705851/2251799813685248

If it were doing an "exact" conversion then it would have ended up with symbolic 3^(1/3), same as

>> sym(3)^(1/3)
ans =
3^(1/3)

sym(1.37) should not be confused with sym('1.37') or vpa(1.37):

>> digits(2)
>> vpa(1.37)
ans =
1.4
>> sym(1.37)
ans =
137/100
>> sym('1.37')
ans =
1.4

vpa() of a numeric expression is equivalent to vpa() of sym() of the numeric expression -- so it does the implicit conversion to extended rational and then calculates the decimal version to as many digits as the current "digits" setting.

sym() of a quoted numeric, on the other hand, treats the quoted value as an exact decimal specification which is not to be approximated:

>> sqrt(2)
ans =
           1.4142135623731
>> sym(1.4142135623731)
ans =
2^(1/2)
>> sym('1.4142135623731')
ans =
1.4142135623731
>> vpa(1.4142135623731)
ans =
1.4142135623730950488016887242097

Now, when I talk about "decimal version" and "exact decimal specification", I mislead a bit. It turns out that the symbolic engine does not use decimal representation of floating point numbers internally -- not even representing as a (integer, log10(denominator)) pair.

I had long long thought that the symbolic engine did its float calculations in decimal, but some experiments I did fairly recently showed me that I was wrong -- tests along the line of finding out that sym('0.1')*10 ~= 1 and the boundary conditions on the accuracy are most easily explained if the symbolic engine is operating in chunks of binary with the number of chunks varying according to the decimal digits (I worked out the limiting ratio but I do not happen to recall what it is.)