vpa(4503599627370491.5)

8 views (last 30 days)
Jon Dattorro
Jon Dattorro on 12 Feb 2018
Edited: Pavel Holoborodko on 5 Oct 2023
vpa(4503599627370491.5) produces 4503599627370496.0 in 2017b. Why?
Further, sym(4503599627370491.5)-sym(4503599627370496) produces 0.

Answers (4)

Walter Roberson
Walter Roberson on 12 Feb 2018
sym and vpa are normal functions. Their arguments are evaluated with the normal MATLAB rules before the function is called. Double precision losses information on a number that large before it gets passed on.
You should enclose large numbers in quotes for use with vpa or syms.
Note: the use of sym on character vectors representing numeric constants will continue to be supported, as will sym on a character string representing a simple variable name. The functionality being discontinued within a small number of releases is using sym() on a character string representing an expression. sym() is defined as treating the character strings as MuPad expressions and the replacement functionality str2sym will instead treat it as MATLAB syntax.

Jon Dattorro
Jon Dattorro on 7 Sep 2023
Advanpix Multiprecision Computing Toolbox (MCT) is far superior to Matlab Variable Precision Arithmetic (VPA) which produces incorrect answers and erroneous results in many different circumstances. This holds true from my own experience and numerical experiments. There is some discussion of that in their Forum here:
I wish MCT were integrated into Matlab because there would be significant improvements to its speed of execution by doing so. It is certainly faster than VPA which requires the Symbolic Math Toolbox.
Mathworks needs to ignore their internal NIH pride and embrace the groundbreaking achievement of Pavel Holoborodko who first introduced MCT in 2011.
Jon Dattorro
  11 Comments
Paul
Paul on 25 Sep 2023
Edited: Paul on 25 Sep 2023
When will vpa return a symfun? Something like this?
syms x
f(x) = x/sym(pi)
f(x) = 
g(x) = vpa(f(x))
g(x) = 
Are there any other cases where vpa returns a symfun?
Walter Roberson
Walter Roberson on 25 Sep 2023
vpa() will return a symfun if the parameter to vpa() evaluates to a symfun .
syms x
f(x) = x/sym(pi)
f(x) = 
g1 = vpa(f)
g1(x) = 
g2 = vpa(3*x + tan(x)^2 + f)
g2(x) = 
g3 = vpa(f^3)
g3(x) = 
g4 = vpa( symfun(gamma(x)^2, x) )
g4(x) = 
g5 = vpa(f(x)) %f(x) is not a symfun
g5 = 
whos
Name Size Bytes Class Attributes cmdout 1x33 66 char f 1x1 8 symfun g1 1x1 8 symfun g2 1x1 8 symfun g3 1x1 8 symfun g4 1x1 8 symfun g5 1x1 8 sym x 1x1 8 sym
vpa() will not return a symfun in any case where the expression does not evaluate to a symfun. In the example of f(x): invoking a symfun() with a definite parameter evaluates the function on the given parameter, returning a symbolic expression.

Sign in to comment.


Pavel Holoborodko
Pavel Holoborodko on 25 Sep 2023
Edited: Pavel Holoborodko on 25 Sep 2023
@"Expecting vpa() to return anything other than sym (or sometimes, symfun) is a substantial misunderstanding of the Symbolic Toolbox...."
I see, the "vpa" is actually a factory function, generating the instances of "sym" class. I stand corrected, "sym.diff" is the full name of the function.
However, these technical details of SMT implementation have nothing to do with the point we are discussing. We discuss the semantics of the functions in SMT. In other words, what the SMT functions do (or should do) with the input arguments. The type they return is just a technicality in this context.
The SMT does much more than pure symbolic operations. Many functions in SMT apply numerical algorithms (using extended precision) to compute the output results. For example, take a look on EIG:
% LAMBDA = EIG(VPA(A)) and [V,D] = EIG(VPA(A)) compute numeric eigenvalues
% and eigenvectors using variable precision arithmetic. If A does not
% have a full set of eigenvectors, the columns of V will not be linearly
% independent.
Same for SVD, or other cases where result cannot be computed analytically or arbitrary-precision plays an important role.
(Of course the functions still return result of type "sym", but the result itself was computed numerically. I hope this distinction is clear.)
This is exactly the case with "diff". Extended precision is important for computing numerical derivatives/differences accurately, since it is a highly ill-conditioned problem (e.g. same as eigenvalues). Therefore it is quite natural to expect the "sym.diff" to have this functionality.
Besides, the original meaning of "diff" stands the same since the MATLAB inception - computing differences and numeric derivatives. It is not clear why SMT developers decided to change its semantics in SMT. They could've just introduced the new function for symbolic derivatives, as @Paul said. Ambiguity would be avoided.
Here is simple example for illustration:
>> A = vpa(rand(3)); % Store numeric matrix as 'sym'
>> ev = eig(A) % Computes e-vals NUMERICALLY using extended precision (see sym.eig)
ev =
1.752670342542888271728211766393
0.8398632587280311835223044391663
-0.18794383321803341008176836388774
>> class(ev) % Result stored as 'sym', but it is inherently numeric
ans =
'sym'
Similarly, I expect the "sym.diff" to provide the analogous functionality of built-in "diff" and to work as "sym.eig" above - compute result numerically using extended precision. But instead we see:
>> d = diff(A)
d =
[0, 0, 0]
[0, 0, 0]
[0, 0, 0]
  11 Comments
Walter Roberson
Walter Roberson on 29 Sep 2023
I do not know when it was introduced. I find it documented at least as far back as R2018b, https://www.mathworks.com/help/releases/R2018b/matlab/matlab_oop/class-support-for-array-creation-functions.html (which is as far as the archives go at the moment)
Pavel Holoborodko
Pavel Holoborodko on 3 Oct 2023
Edited: Pavel Holoborodko on 5 Oct 2023
This covers only the simplest case when the user passes the typename argument explicitly (e.g. zeros(...,'mp')).
However, in most of the cases, existing/legacy code relies on default form, without the typename argument.
That is why our users requested to provide a special mode in which the array-creation functions generate extended-precision arrays by default (e.g. when no typename is specified, or it refers to standard floating-point types - 'double' or 'single').
This is useful when legacy code of considerable volume needs to be converted to extended precision:
>> mp.OverrideDoubleBasicArrays(true); % Basic arrays are created as 'mp' by default
>> A = zeros(3);
>> whos A
Name Size Bytes Class Attributes
A 3x3 272 mp
>> % Run legacy code with zeros(...), eye(...), etc.
>> % ...
>> mp.OverrideDoubleBasicArrays(false); % Return to default behavior
That is why our toolbox has to supersede the default functions and call them only when applicable (using 'builtin' calls). This approach works well, with the exception of slow 'builtin' calls.
***
Another interesting case is when the functions are called without arguments at all, e.g. zeros, eye, etc. This case is treated differently by MATLAB - it just returns 0 or 1, without even calling the global/built-in functions. No OOP involved, this is just hardcoded in MATLAB. Meaning these cases are impossible to overload for custom class types at all. Even by superseding the built-in functions.

Sign in to comment.


Jon Dattorro
Jon Dattorro on 3 Oct 2023
Respectfully:
Are there any further objections to tight embedding of Advanpix Multiprecision Computing Toolbox into core Matlab functionality?
Pavel has made one concession (in documentation) and is willing to make another (in coding per Mr Lord).
Jon Dattorro

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!