Multiply two variables in Matlab with vpa - high precision

Question

I want to be sure that two variables, a and b, are multiplied with high precision, i.e., perform the product c of a and b with arbitrary precision (in my case 50 correct decimal digits):

a = vpa(10/3,50);  
b = vpa(7/13,50);  
c = eval(vpa(vpa(a,50)*vpa(b,50),50));      % I basically want to do just c = a*b;    

which gives me

a = 3.3333333333333333333333333333333333333333333333333
b = 0.53846153846153846153846153846153846153846153846154
c = 23.333333333333333333333333333333

Testing

d = eval(vpa(c*13,50))

gives

d = 23.333333333333333333333333333333333333335292490585

which shows that the multiplication to get c was not carried out with 50 significant digits.

What's wrong here, but, more importantly, how do I get a correct result for a*b and for other operations such as exp?

Answer 1

First, should use vpa('7/13',50) or vpa(7,50)/13 to avoid the possibility of losing precision dues to 7/13 being calculated in double precision floating point (I believe that vpa, like sym, tries to guess common constants and rational fractions, but you shouldn't rely on it).

The issue is that while a and b are stored as 50-digit variable precision values, your multiplication is still being performed according to the default value of digits (32). The second argument to vpa only appears to specify the precision of the variable, not any subsequent operations on or with it (the documentation is not particularly helpful in this respect).

One way to accomplish what you want would be:

old_digits = digits(50);
a = vpa('10/3')
b = vpa('7/13')
c = a*b
d = 13*c
digits(old_digits);

Another would be to use exact symbolic expressions for all of the math (potentially more expensive) and then convert the result to 50-digit variable precision at the end:

a = sym('10/3')
b = sym('7/13')
c = a*b
d = vpa(13*c,50)

Both methods return 23.333333333333333333333333333333333333333333333333 for d.