Precise conversions from double to symbolic

I have the following number

r = 1.78503

I want to obtain the exact symbolic representation of 1/r, and I applied

p = 1/sym(r);

which gives me 1125899906842624/2009765110711289.

However, if I apply the form sym(1/r) then I got 2522982598259131/4503599627370496 which is different from the previous one. I understand this is due to the floating point numeric 1/r so the later form may not be accurate. Based on this observation, what measures should be taken in order to have exact values in a situation like this ? Is the form 1/sym(r) fully able to extract the exact symbolic representation here ? Thanks.

[James Tursa]

Can you clarify this post?

In the first place, 2.1234 cannot be represented exactly in IEEE double:

>> r = 2.1234
r =
    2.1234
>> num2strexact(r)
ans =
2.123400000000000176214598468504846096038818359375

So when you do this:

>> r = 2.1234
r =
    2.1234
>> sym(r)
ans =
10617/5000

MATLAB is actually making some assumptions behind the scenes that you really meant the exact decimal value 2.1234 to be converted and not the actual number that is present in the IEEE double value.

E.g., the trailing bits have to be different from the closest representation of 2.1234 by 7 bits before this assumption is tossed and an exact conversion of the exact bit pattern is calculated:

>> sym(r+2^6*eps(r))
ans =
10617/5000
>> sym(r+2^7*eps(r))
ans =
597683965547423/281474976710656

And, as you can see, I am getting different numbers than you are showing above for the conversions. The numbers you are showing give:

>> 1125899906842624/2009765110711289
ans =
    0.5602
>> 1/ans
ans =
    1.7850
>> 2522982598259131/4503599627370496
ans =
    0.5602
>> 1/ans
ans =
    1.7850

So, what are you really working with here? And why do you need "exact" conversions? What do you mean by "exact"? What are you doing with this downstream?

[Qian Feng]

My apology for the mistake in my code, the value of r should be 1.78503.

[Walter Roberson]

Change the assignment in the code I showed, and run the steps, and pick the version that you want. 2522982598259131/4503599627370496 or 1125899906842624/2009765110711289 or 100000/178503 or 560214674263177649675355596264489/1000000000000000000000000000000000 (the last would change if you changed the number of digits you have set.)

Answer 1

[Walter Roberson]

Compare:

r = 2.1234
sym(1/r)
1/sym(r)
1/sym(r,'d')
1/sym(r,'e')
1/sym(r,'f')
1/sym(r,'r')

s = sym(sprintf('%.16g', r));
feval(symengine,'numeric::rationalize',1/s,'Exact')
1/feval(symengine,'numeric::rationalize',s,'Exact')

[Qian Feng]

BTW, is there a possible way to convert a vpa object into symbolic rational representation ? Using sym() does not work.

[Walter Roberson]

feval(symengine, 'numeric::rationalize', TheVpaObject, 'Exact')

[Karan Gill]

Hold on! You should not be converting vpa into exact representation. In your workflow, you are supposed to only go from symbolic to numeric forms. Why are you going the other way?

[Walter Roberson]

Oh, I do this all the time. When people post asking for solutions to equations, the majority of the time they post using floating point numbers, and to find the exact solution I have to rewrite the floating point numbers as rationals to put it through (a different vendor's) symbolic routines. Even people who are not doing anything wonky like raising a value to a floating point power tend to write their expressions with constants like 0.25 or 0.01 or 0.3333, but more common is that arbitrary values like 8.6047E-26 will be in the expressions. It is a continued struggle to point out to people that they are trying to make calculations that cannot be mathematically justified because of the low precision of their inputs -- cases where, for example, a change between 8.60471E-26 8.60473E-26 might completely change the sign of a value that they claim that they need the "exact" solution for.

With the Symbolic Toolbox, if you mix symbolic and floating point then by default the floating point is converted as if by sym() with the 'r' flag, but that is documented as being an approximate conversion, so you really need to know what the alternatives are, in case you want the 8.6047E-26 to convert to sym('86047/10^30') equivalent (so that the number used reflects exactly what it looks like visually... if only because you want to proceed to prove that using the exact value the user specified leads to numeric nonsense.)

[Karan Gill]

If you want to pass a number "exactly" then using string input to sym is okay: https://www.mathworks.com/help/symbolic/sym.html#bu7u7ur-6

>> sym('1.78503')
ans =
1.78503
>> vpa(sym('1.78503'))
ans =
1.78503

Does this help?

[Qian Feng]

Thanks Karan for the suggestion. My question here is related to another question I asked in https://au.mathworks.com/matlabcentral/answers/311003-numerical-values-of-an-non-elementary-integral. For this case, it is better to have vpa(E1) first instead of double(vpa(E1)) since E1 still involves some calculations subsequently. For the aforementioned case I am dealing with vpa is sufficient in fact. However, I would like to know the means to transfer vpa into rational representation which might be applied to other potential problems, since it will not suffer any round error. My question here is what happen when an operation involves in both rational and vpa representations ? Does it transfer the rational into vpa automatically ? If so, that is in fact the reason why I asked about converting vpa into exact representations.

[Walter Roberson]

When you mix floating point and symbolic, the floating point is converted as if by sym() with the 'r' flag, which is not an exact transformation. For example since 1/3 is not exactly representable in binary floating point, for sym(1/3) to produce sym(1)/sym(3) requires approximation of the 0.333333333333333314829616256247390992939472198486328125 (exact value of 1/3 in floating point) .

I listed the possible conversion options above, the 'd', 'e', 'f', and 'r' flags.

[Qian Feng]

Yes I understand that now. To get the exact value of 1/3 we can sym(1)/sym(3) and it does not have a error problem in this case. How about if we mix 'r' representation with 'd' objects ? I found out it seems that the rational object is automatically transferred into digits form.

[Walter Roberson]

Correct, when there are rational and decimal numeric constants being added or multiplied with each other, the rational are converted to decimal. Also, elementary functions such as exp() and cos() will be evaluated when a decimal format sym is passed in, but not when a rational is passed in. (powers of a rational might be simplified but not fully evaluated.)

Answer 2

[John D'Errico]

Too late of course. But the point is that when you create r as a double:

r = 1.78503;

Then it is NOT stored exactly, as 1.78503. Nothing you do will then allow MATLAB to know that you really intended 1.78503, and not the number actually stored, which is...

sprintf('%0.55f',1.78503)
ans =
1.7850299999999998945554580132011324167251586914062500000

Even if you try to pass that number into sym, MATLAB will get it wrong, because you passed in a double precision number as far as sym was concerned.

vpa(sym(1.78503),55)
ans =
1.78502999999999989455545801320113241672515869140625

A simple solution is to go directly to symbolic form, but even there one must be careful.

r = sym('1.78503')
r =
1.78503
vpa(r,55)
ans =
1.785029999999999999999999999999999999999842248658119649

So it looks like r only had about 40 decimal digits stored. (As a guess, roughly 128 bits in the mantissa.) Still better than 16 digits though.

You can do better, by avoiding decimals completely. Integers work best.

r = sym('178503/100000')
r =
178503/100000
vpa(r,300)
ans =
1.78503

Or, you can use my HPF toolbox.

hpf('1.78503',100)
ans =
1.78503

[James Tursa]

OP still hasn't stated why there is a need for "exact" conversions and what they are being used for downstream. So I have yet to be convinced that this all isn't just a pointless exercise ...

Answer 3

[Karan Gill]

Use quotes to keep your input exactly as it.

r = sym('1.78503')

But you don't get the fractional representation since you asked to keep it exactly as the decimal.

>> p = 1/r
p =
0.56021467426317764967535559626449