4.1370454708904844665084965527057647705078125
Why is vpa not converting double accurately?
Show older comments
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)
sf = sprintf('%.50f',f)
str2sym(sf)
vpa(f)
double(vpa(f)) - f
vpa(f,50) % Using this form of function call doesn't help
digits 200 % this doesn't help either
vpa(f)
double(vpa(f)) - f % back conversion to double doesn't match either
What is going on here? Obviously vpa can store the exact floating point binary to decimal conversion of f:
vpa(str2sym(sf))
double(ans) - f
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
fprintf('%.50f\n',r)
vpa(r)
r = rand*100-50
fprintf('%.50f\n',r)
vpa(r)
r = rand*100-50
fprintf('%.50f\n',r)
vpa(r)
Accepted Answer
John D'Errico
on 22 Mar 2025
Edited: John D'Errico
on 22 Mar 2025
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)
num2hex(f)
Good. That works.
vf = vpa(f)
num2hex(double(vf))
and clearly that is not the same number. So, is the problem in VPA? Instead, I'll try this.
sym(f)
Lol. That is sort of interesting. And a bit of a surprise.
digits 200
vpa(sym(f))
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. ;-)
19 Comments
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)
digits
double(vpa(sym(f,'d'))) - f
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'))
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.
format long
digits 50
hf = '40108c55a5de2880'; % making sure I am starting from the value I think I am
f = hex2num(hf)
sym(f)
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])
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 )
sym( sqrt(241) * sqrt(16499) ) / 482
sym( sqrt(241 * 16499) ) / 482
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)
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'))
Walter Roberson
on 22 Mar 2025
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
on 22 Mar 2025
Edited: James Tursa
on 22 Mar 2025
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?
format long g
digits 25
f = 1/3
fprintf('%.999g\n', f)
sym(f,'d')
sym(f,'e')
vpa(ans)
subs(ans, sym('eps'), eps)
sym(f,'f')
vpa(ans)
sym(f,'r')
vpa(ans)
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
on 22 Mar 2025
Edited: John D'Errico
on 22 Mar 2025
"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
on 22 Mar 2025
Edited: James Tursa
on 22 Mar 2025
Paul
on 22 Mar 2025
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))
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);
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
fprintf('%d discrepenciesB found\n', counter);
John D'Errico
on 22 Mar 2025
Interesting there were only a few discrepancies.
Walter Roberson
on 23 Mar 2025
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.
format long
digits 50
hf = '40108c55a5de2880'; % making sure I am starting from the value I think I am
f = hex2num(hf)
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')
should be exactly equal in value to f. But it isn't (or at least appears not to be)
avpa = vpa(asym)
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
on 23 Mar 2025
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.
"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)
which matches the SYM value exactly. I believe that up to 1075 digits are required, in the worst case.
format long g
E = eps(0)
S = sprintf('%.2000f', E);
S = regexprep(S, '0+$', '')
length(S)
The 1076 includes the leading 0 and the decimal place
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
digits 20
vpa(0.1) + vpa(0.2) - vpa(0.3)
digits 21
vpa(0.1) + vpa(0.2) - vpa(0.3)
digits 22
vpa(0.1) + vpa(0.2) - vpa(0.3)
digits 23
vpa(0.1) + vpa(0.2) - vpa(0.3)
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)
2^(-100)
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)
digits 100
vpa(ans)
fprintf('%.100f\n',2^(-94))
And yes, you can see that the residual value stored is in fact exactly the value of a binary bit.
Walter Roberson
on 13 May 2025
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.
More Answers (0)
Categories
Find more on Code Performance in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
