Why is vpa not converting double accurately?

Question

James Tursa 2025-3-22

3
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2175497-why-is-vpa-not-converting-double-accurately

编辑： James Tursa 2025-5-15

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

John D'Errico 2025-3-22

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2175497-why-is-vpa-not-converting-double-accurately#answer_1562234

编辑：John D'Errico 2025-3-22

在 MATLAB Online 中打开

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. ;-)

19 个评论
显示 17更早的评论隐藏 17更早的评论

John D'Errico 2025-3-22

在 MATLAB Online 中打开

@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 2025-3-22

编辑：James Tursa 2025-3-22

在 MATLAB Online 中打开

@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

James Tursa 2025-3-22

编辑：James Tursa 2025-3-22

在 MATLAB Online 中打开

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

Walter Roberson 2025-3-22

在 MATLAB Online 中打开

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

Walter Roberson 2025-3-23

在 MATLAB Online 中打开

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 2025-5-13

编辑：James Tursa 2025-5-15

在 MATLAB Online 中打开

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 2025-5-13

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

请先登录，再进行评论。

Why is vpa not converting double accurately?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

19 个评论
显示 17更早的评论隐藏 17更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

Why is vpa not converting double accurately?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

19 个评论 显示 17更早的评论隐藏 17更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

19 个评论
显示 17更早的评论隐藏 17更早的评论