Why is vpa not converting double accurately?

I ran across this example that I don't understand. I don't use the Symbolic Toolbox much, but I thought vpa by default converted doubles to symbolic representations that were an accurate conversion up to so many decimal digits. I.e., an exact floating point binary to decimal conversion up to the specified number of decimal digits. But the following did not match my expectations:

format long
digits 50
hf = '40108c55a5de2880'; % making sure I am starting from the value I think I am
f = hex2num(hf)

f = 
   4.137045470890484

sf = sprintf('%.50f',f)

sf = '4.13704547089048446650849655270576477050781250000000'

str2sym(sf)

ans = 

4.1370454708904844665084965527057647705078125

vpa(f)

ans = 

4.1370454708905203952304813814278565783465431783307

double(vpa(f)) - f

ans = 
     3.552713678800501e-14

vpa(f,50) % Using this form of function call doesn't help

ans = 

4.1370454708905203952304813814278565783465431783307

digits 200 % this doesn't help either
vpa(f)

ans = 

4.1370454708905203952304813814278565783465431783307106731417481082606305993533207364559999183814049903246844172718745791741538190785178147855046143279189519799161307279637872963853860781001401307209969

double(vpa(f)) - f % back conversion to double doesn't match either

ans = 
     3.552713678800501e-14

What is going on here? Obviously vpa can store the exact floating point binary to decimal conversion of f:

vpa(str2sym(sf))

ans = 

4.1370454708904844665084965527057647705078125

double(ans) - f

ans = 
     0

I have to use fprintf with str2sym to get an accurate conversion. Why doesn't vpa do that with f in the first place? Is there some setting I am unaware of? And what are all those digits pouring out of vpa for this f?

Is this some type of p/q morphing going on in the background?

Other random values generally match my expectations. E.g.,

r = rand*100-50

r = 
  36.006640442455534

fprintf('%.50f\n',r)

36.00664044245553441214724443852901458740234375000000

vpa(r)

ans = 

36.00664044245553441214724443852901458740234375

r = rand*100-50

r = 
 -47.047616571398464

fprintf('%.50f\n',r)

-47.04761657139846420250250957906246185302734375000000

vpa(r)

ans = 

r = rand*100-50

r = 
   4.524320317487962

fprintf('%.50f\n',r)

4.52432031748796248393773566931486129760742187500000

vpa(r)

ans = 

4.524320317487962483937735669314861297607421875

Answer 1

[John D'Errico]

Quite interesting. When I first saw the title, I assumed it would be an obvious mistake. Then I saw the poster, and knew it would be interesting. My first assumption is the problem must lie in VPA, which sometimes seems to stand for:

VPA - Very Poor Arithmetic

If I look at what is happening though...

format long
digits 50
hf = '40108c55a5de2880'; % making sure I am starting from the value I think I am
f = hex2num(hf)

f = 
   4.137045470890484

num2hex(f)

ans = '40108c55a5de2880'

Good. That works.

vf = vpa(f)

vf = 

4.1370454708905203952304813814278565783465431783307

num2hex(double(vf))

ans = '40108c55a5de28a8'

and clearly that is not the same number. So, is the problem in VPA? Instead, I'll try this.

sym(f)

ans = 

Lol. That is sort of interesting. And a bit of a surprise.

digits 200
vpa(sym(f))

ans = 

4.1370454708905203952304813814278565783465431783307106731417481082606305993533207364559999183814049903246844172718745791741538190785178147855046143279189519799161307279637872963853860781001401307209969

Hmm. Is that really what is happening? What does that strange number with the radicals resolve to as a decimal? I'll try HPF, since I know exactly what HPF is doing. Hey, I trust HPF implicitly, since I know the guy who wrote it.

DefaultNumberOfDigits 200
x = sqrt(hpf(241))*sqrt(hpf(16499))/hpf(482)
x =
4.1370454708905203952304813814278565783465431783307106731417481082606305993533207364559999183814049903246844172718745791741538190785178147855046143279189519799161307279637872963853860781001401307209969

Ah. So now I know what happened. Well, I think I do, I think I do.

When you do this:

vf = vpa(f)

MATLAB FIRST converts the number to a sym, because VPA only understands symbolic numbers. Effectively, it expands it as

vf = vpa(sym(f))

But what does sym do? It decides the result of sym(f) is sqrt(241)*sqrt(16499)/482. After all, that is what we would all do, right? The obvious form. Then, and ONLY then, does it throw it into VPA.

So the problem is not in VPA. The problem lies in the decision to approximate that floating point number as a sym in the form it chose. The problem lies in sym.

My head hurts, just a little. ;-)

[Paul]

I'm not so sure the problem lies with sym, though I do think that there is a gap at vpa and a potential inconsistency between the doc pages sym and vpa wrt to how those functions handle numeric inputs.

Assuming that vpa does call sym (or its equivalent) under the hood, it appears that it's using the default flag input to sym, which is the 'r' option, which "might not represent the floating-point value precisely."

format long
digits 50
hf = '40108c55a5de2880'; % making sure I am starting from the value I think I am
f = hex2num(hf);
a1 = double(vpa(f)) - f;           % whatever it does
a2 = double(vpa(sym(f))) - f;      % explicit conversion with default flag
a3 = double(vpa(sym(f,'r'))) - f;  % explicit conversion with explicit flag
isequal(a1,a2,a3)

ans = logical
   1

The 'd' flag for sym works as expected with the default digits

digits

 
Digits = 50
 

double(vpa(sym(f,'d'))) - f

ans = 
     0

[John D'Errico]

@Paul

I'd argue it is an issue with sym. Why? Because the result that VPA returns is EXACTLY:

sqrt(241)*sqrt(16499)/482

I showed that when I computed that same number using HPF, to as many digits in both cases. And that means VPA is seeing the number with the radicals, and not the original HEX number as a floating point value.

If you want to argue the issue lies in how sym is called by VPA, you can. But the point is, VPA is not seeing the floating point number, but the radical-ised approximation we saw sym return. That is the fundamental problem. VPA did what it is designed to do.

And, yes, you can force sym to use the number itself. Looks like you show that yourself.

format long
digits 50
hf = '40108c55a5de2880'; % making sure I am starting from the value I think I am
f = hex2num(hf);
vpa(sym(f,'d'))

ans = 

4.1370454708904844665084965527057647705078125

What does this mean to me? I assume vpa calls sym with the default, but the default when VPA sees a floating point number should NOT be to convert the number to an approximation FIRST. So if you want to say this is not an issue with sym but perhaps in how sym was called, feel free, though it seems a silly distinction to me.

Are you happier if we decide the problem lies in the interface between the two? So is the problem that SYM is making the wrong choice by default, or should it be that VPA should force SYM to make a different choice?

Should SYM be smarter? Hard to say. Should VPA be smarter? This too is hard to say. For example, suppose we were to change VPA, so that when we call

vpa(FLOAT)

where float is some double precision number, then it automatically uses the 'd' flag? Now @James Tursa, or someone else, possibly you, might have noticed that

vpa(sym(FLOAT))

returns a different result. Surely that too would be a reportable inconsistency.

What does this mean to me then? The problem lies in SYM! This is because if you would want vpa(FLOAT) to encode the exact floating point number, then surely vpa(sym(FLOAT)) should return the same number.

[James Tursa]

@John D'Errico

format long
digits 50
hf = '40108c55a5de2880'; % making sure I am starting from the value I think I am
f = hex2num(hf)

f = 
   4.137045470890484

sym(f)

ans = 

OH! Duh! Of course, why didn't I think of this? (slaps hand on my head). This is obviously one of those "Common Double-Precision Inputs" mentioned in the vpa doc! I feel embarrassed and a bit silly for even posting this question in the first place. I'm sure Ramanujan would have recognized this immediately, and John D'Errico shortly thereafter, but sadly I am not of that caliber. (And my head hurts a little.)

On a more serious note, I am glad to learn what is going on. In particular, the fact that there is slop in the background that creates "pockets" of double precision numbers that convert to the same sym value. E.g.,

sym(f+eps(f)*[-12;-11;0;91;92])

ans = 

This is all, I presume, to help the naive programmer. I mean, I get it. If I worked at The Mathworks, would I want to be fielding thousands of complaints from new programmers about sym or vpa being broken for returning something like this:

0.333333333333333314829616256247390992939472198486328125

for a simple input like 1/3? Of course not. That is the understandable motivation. But it doesn't really help the seasoned programmer at all. If I caught myself coding something like sym(1/3) instead of sym(1)/sym(3), I would rap my knuckles with a ruler to ensure I didn't do it again. It would be poor programming practice in my book. I would never rely on the sym engine to fix my bad programming. More importantly, I think it renders direct vpa(floating_point) expressions practically unusable for any precision work, because the results are unpredictable for arbitrary inputs. You get stuff like this happening because the sym engine obviously has limits for what it can recognize:

sym( sqrt(241) * sqrt(16499) / 482 )

ans = 

sym( sqrt(241) * sqrt(16499) ) / 482

ans = 

●

sym( sqrt(241 * 16499) ) / 482

ans = 

Three different answers depending on where the arithmetic takes place. Again, I would rap my knuckles with a ruler if I coded any of the above. The proper way to code it would be to convert the integers to sym first. E.g., something like this would suffice:

sqrt(sym(241)) * sqrt(sym(16499)) / sym(482)

ans = 

Bottom line for me is to never use vpa(floating_point) because I can never be sure what it will convert to in the background. I guess I will always use vpa(sym(floating_point,'d')) going forward if I want to use vpa. E.g., this is what I would typically want as a result when doing extended precision work with vpa:

vpa(sym(f,'d'))

ans = 

4.1370454708904844665084965527057647705078125

[Walter Roberson]

I would argue that you should be using

vpa(sym(f,'f'))

The f flag signals that the exact rational equivalent of the double precision number should be created.

[James Tursa]

@Walter Roberson Noted, but is there a case where the vpa result would be different? Is vpa(sym(f,'d')) redundant? I.e., do I really even need the outer vpa( ) call for the 'd' option? Put another way, is sym(f,'d') equivalent to vpa(sym(f,'f'))?

[Paul]

@James Tursa

FWIW, I don't think that outer vpa is needed when using the 'd' option and I should not have included it in this comment. Oversight on my part. Would be very interested if there is a counterexample.

The doc page for vpa agrees with you: "For accurate results, convert numeric expressions to symbolic expressions with sym [on input to vpa]." Of course, the user is then obligated to use the appropriate flag input to sym and/or do the sym arithmetic appropriately, as you showed above.

Do you think that either vpa or sym demonstrate actual behavior that is counter to the documentation?

[Walter Roberson]

format long g
digits 25
f = 1/3

f = 
         0.333333333333333

fprintf('%.999g\n', f)

0.333333333333333314829616256247390992939472198486328125

sym(f,'d')

ans = 

0.3333333333333333148296163

sym(f,'e')

ans = 

vpa(ans)

ans = 

subs(ans, sym('eps'), eps)

ans = 

0.3333333333333333148296163

sym(f,'f')

ans = 

●

vpa(ans)

ans = 

0.3333333333333333148296163

sym(f,'r')

ans = 

vpa(ans)

ans = 

0.3333333333333333333333333

So at least in this test, sym(f,'d') works out the same as vpa(sym(f,'f')).

But sometimes you want f = 1/3 (double precision) to be converted as 1/3 rational, and sometimes you want f = 1/3 (double precision) to be converted to the exact rational fraction ...

[John D'Errico]

"I'm sure Ramanujan would have recognized this immediately, "\

Probably so. But I wonder how sym chose this approximation, as it is not at all obvious why that one would have been chosen over a milllion others, probably equally as good. Yes, my head hurts :)

Anyway, I think the problem is, that sometimes you want sym (and vpa) to make different choices. Good to have the flag in sym be available. But I would prefer that vpa(float) ALWAYS return by default the exact floating point number.

[James Tursa]

@Paul I wouldn't necessarily say that the vpa or sym doc is wrong when it comes to vpa(floating_point), but it is certainly lacking in clarity. For example, in the vpa description the only allowed input mentioned up front is that x must by symbolic. So vpa(floating_point) isn't even mentioned here.

Scanning through the doc the examples are at first careful to convert floating point to sym first before calling vpa. E.g., this

p = sym(pi);
pVpa = vpa(p)

instead of this

pVpa = vpa(pi)

But later on the doc starts to use direct floating point input, e.g.,

pVpa = vpa(pi,100)

Here there is a double precision input, and there is an implicit assumption about what is going on in the background, but the doc doesn't point out this distinction and the sudden change in input type.

It's not until you get to the section "Restores Precision of Common Double-Precision Inputs" that you get the clues for how the previous examples involving floating point input in the doc worked. The fact that there are pockets 100*eps wide (at least) that all convert to the same value (part of the "round-off errors" I presume) isn't stressed at all, while I consider this to be an extremely important fact for how direct floating point conversions actually work.

And frankly, using examples converting 1/3 and 1/10 instead of sym(1)/sym(3) and sym(1)/sym(10) is a disservice to the reader. The fact that the doc is actually encouraging the reader to program this way is bad advice IMO, and completely ignores the pitfalls that can be encountered when doing floating point arithmetic prior to sym conversion.

[Paul]

100% agree.

The description really should be updated to include numeric input.

The term "restores precision" is troublesome. Seems like it means "restores precision to what we think you meant."

The inconsistency between vpa and sym that I reference above is that the former "restores precision" for:

"numeric inputs that match the forms p/q, pπ/q, (p/q)^(1/2), 2^q, and 10^q, where p and q are modest-sized integers."

whereas the default behavior for the latter:

"converts floating-point numbers obtained by evaluating expressions of the form p/q, p*pi/q, sqrt(p), 2^q, and 10^q (for modest-sized integers p and q) to the corresponding symbolic form."

The third elements in those respective lists are different, which suggests a possbility that there could be a case where vpa(f) does not equal vpa(sym(f))

[Walter Roberson]

sym(f, 'd') turns out to be different than vpa(sym(f,'f')), in the digits beyond what is displayed. To the digits displayed the results are identical.

format long g
digits 25
rng(1, 'swb2712')

counter = 0;
for K = 1 : 1000
    f = rand();
    t1 = sym(f, 'd');
    t2 = vpa(sym(f, 'f'));
    ch1 = char(t1);
    ch2 = char(t2);
    if ~isequal(ch1, ch2)
        fprintf('round %d discrepency\n', K);
        fprintf('%s\n', ch1);
        fprintf('%s\n', ch2);
        counter = counter + 1;
    end
end
fprintf('%d discrepencies found\n', counter);

0 discrepencies found

counter = 0;
for K = 1 : 3
    f = rand();
    t1 = sym(f, 'd');
    t2 = vpa(sym(f, 'f'));
    ch1 = char(t1);
    ch2 = char(t2);
    if ~isequal(t1, t2)
        fprintf('\nround %dB discrepency\n', K);
        fprintf('%s\n', ch1);
        fprintf('%s\n', ch2);
        fprintf('difference: %s\n', char(t1-t2));
        counter = counter + 1;
    end
end

round 1B discrepency

0.5407167691652339236441094

0.5407167691652339236441094

difference: 1.828795820512754948543946e-26

round 2B discrepency

0.4082336336977269164449922

0.4082336336977269164449922

difference: 1.544833443857338942762579e-26

round 3B discrepency

0.2083801804355891185149829

0.2083801804355891185149829

difference: -4.254270149539512092510197e-26

fprintf('%d discrepenciesB found\n', counter);

3 discrepenciesB found

[John D'Errico]

Interesting there were only a few discrepancies.

[Walter Roberson]

@John D'Errico

My second for loop only goes to 3. In my tests with larger loops, every pair of values was different. I didn't think it was interesting to display that many cases, but at the same time I thought it was worth running the first loop numerous times to illustrate that the printed values are always the same.

[Paul]

Returning to @James Tursa's example

format long
digits 50
hf = '40108c55a5de2880'; % making sure I am starting from the value I think I am
f = hex2num(hf)

f = 
   4.137045470890484

According to sym, if I use the 'f' option: "The returned symbolic number is a precise rational number that is equal to the floating-point value."

I thought that in this case the "floating-point value" of f is exactly 4.137045470890484 because of way that f was constructed. Therefore sym(f,'f')

asym = sym(f,'f')

asym = 

should be exactly equal in value to f. But it isn't (or at least appears not to be)

avpa = vpa(asym)

avpa = 

4.1370454708904844665084965527057647705078125

What is the source of all of those trailing digits after the 0484? Am I misunderstanding the quoted passage from the doc page (or some other aspect of this example)?

[Walter Roberson]

@Paul

The 'f' option essentially decomposes the IEEE 754 representation into sign, exponent, and mantisa, and constructs the explicit representation, roughly (-1)^sign * 2^(exponent - bias) * (1 + mantissa/2^52)

Remember that floating point numbers are stored in binary, not decimal.

[Stephen23]

"Am I misunderstanding the quoted passage from the doc page (or some other aspect of this example)?"

15 digits are not enough to show the exact numeric value. The numeric value is precisely

hf = '40108c55a5de2880'; % making sure I am starting from the value I think I am
f = hex2num(hf);
fprintf('%.999g\n',f)

4.1370454708904844665084965527057647705078125

which matches the SYM value exactly. I believe that up to 1075 digits are required, in the worst case.

[Walter Roberson]

format long g
E = eps(0)

E = 
     4.94065645841247e-324

S = sprintf('%.2000f', E);
S = regexprep(S, '0+$', '')

S = '0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000004940656458412465441765687928682213723650598026143247644255856825006755072702087518652998363616359923797965646954457177309266567103559397963987747960107818781263007131903114045278458171678489821036887186360569987307230500063874091535649843873124733972731696151400317153853980741262385655911710266585566867681870395603106249319452715914924553293054565444011274801297099995419319894090804165633245247571478690147267801593552386115501348035264934720193790268107107491703332226844753335720832431936092382893458368060106011506169809753078342277318329247904982524730776375927247874656084778203734469699533647017972677717585125660551199131504891101451037862738167250955837389733598993664809941164205702637090279242767544565229087538682506419718265533447265625'

length(S)

ans = 
        1076

The 1076 includes the leading 0 and the decimal place

[James Tursa]

There are clues in the examples posted in this thread, which led me to an investigation. I have now convinced myself that vpa( ) in the background stores numbers in an extended binary (or hex?) floating point format. I'm guessing the value itself uses a whole number of bytes or words. So when you set digits or convert a number using a certain number of requested digits, what actually must happen in the background is a calculation for how many binary bits are needed in a significand to ensure the requested number of decimal digits is met (and add in room for sign and exponent). Then extend that number to get you on a byte or word boundary, and add guard bytes or words. And that is what you actually end up with in the background. All vpa arithmetic in the background is probably binary (or equivalent) floating point arithmetic. And, if that is the case, then we should expect to see binary floating point artifacts in the results. And we do. E.g.,

format longg
0.1 + 0.2 - 0.3

ans = 
      5.55111512312578e-17

digits 20
vpa(0.1) + vpa(0.2) - vpa(0.3)

ans = 

5.0487097934144755546e-29

digits 21
vpa(0.1) + vpa(0.2) - vpa(0.3)

ans = 

0.0

digits 22
vpa(0.1) + vpa(0.2) - vpa(0.3)

ans = 

digits 23
vpa(0.1) + vpa(0.2) - vpa(0.3)

ans = 

0.0

The double precision calculation is a common one to use for this demonstration. But note that vpa has some very similar artifacts, depending on how many decimal digits were requested, which I presume affected the actual number of bytes/words used in the background.

The fact that the residuals are precisely the value of a binary bit lends credence to this hypothesis:

2^(-94)

ans = 
      5.04870979341448e-29

2^(-100)

ans = 
      7.88860905221012e-31

I find it interesting that the number of effective binary bits where the residual occurred in the second case was 6 bits away from the first case instead of 8.

Examining this example in more detail:

digits 20
vpa(0.1) + vpa(0.2) - vpa(0.3)

ans = 

5.0487097934144755546e-29

digits 100
vpa(ans)

ans = 

0.0000000000000000000000000000504870979341447555463506281780983186990852118469774723052978515625

fprintf('%.100f\n',2^(-94))

0.0000000000000000000000000000504870979341447555463506281780983186990852118469774723052978515625000000

And yes, you can see that the residual value stored is in fact exactly the value of a binary bit.

[Walter Roberson]

@James Tursa

If it uses the same strategy as Maple does (MuPAD was designed as a Maple work-alike) then the strategy is something like using chains of 32 bit integers, out of which only 30 bits are (normally) used. Maple does not use base 100,000,000 (a number slightly below 2^30); it uses base 2^30 if I recall correctly.

At the moment, I do not recall how the decimal place information is stored.