I'm having a bit of an issue regarding numerical precision and I'm not sure how to deal with it.
I have a certain randomly generated matrix, say $M$, that I wish to compute the eigenvalues. The entries of this matrix are real and span many orders of magnitude, from $10^{16}$ down to $10^{-14}$.
I am particularly interested in the case where the smallest eigenvalue is bounded in a region, say $10^{-10}<\lambda_{min}<10^{-9}$. My program then generates some random numbers for the entries of $M$, computes the eigenvalues through
EigenValues[N[M,precision]]
where 'precision' is specified somewhere before the above evaluation$^a$, and then stores some of the parameters for $M$ that meets the criteria above.
In order to test the numerical precision of my numbers I generated an example instance of $M$, let's call it $M0$. I write down its entries in exact form, i.e.
Precision[M0]=Infinity
I then compute the eigenvalues with different precisions, namely consider the two cases
Eigenvalues[N[M0,MachinePrecision]] Eigenvalues[N[M0,30]]
and test these results against the same computation carried out by a C program.
For a specific matrix $M0$ that I am using for this test I get the following results for its eigenvalues.
Mathematica With MachinePrecision: $$ 1.08568\times 10^{15} \\ -1.08568\times 10^{15} \\ 2.04979\times 10^{13} \\ 1.98037\times 10^{13} \\ -1.51195\times 10^{10} \\ 1.51195\times 10^{10} \\ 737.092 \\ 709.496 \\ -288.09 \\ 287.764 \\ -0.0527701 $$
Mathematica With precision=30: $$ 1.08568468076658564504709826065\times 10^{15} \\ -1.08568468076544177351940777511\times 10^{15} \\ 2.04978825179076415641037942687\times 10^{13} \\ 1.98037159026622180001959563609\times 10^{13} \\ -1.51195073626019461804045801802\times 10^{10} \\ 1.51195073626019461804034605668\times 10^{10} \\ 737.067283246721413545966471801 \\ 709.473111591071399129829176386 \\ -288.009491702297922120741501773 \\ 287.845096874967906231858249555 \\ -\text{1.04616771734909972710824504576458020590918184213753$\grave{ }$30.*${}^{\wedge}$-8} $$
Custom made C code, where the numbers are taken to be Double $$ +2.472\times10^{-10}\\ +2.878\times10^{+02}\\ -2.880\times10^{+02}\\ +7.095\times10^{+02}\\ +7.371\times10^{+02}\\ +1.512\times10^{+10}\\ -1.512\times10^{+10}\\ +1.980\times10^{+13}\\ +2.050\times10^{+13}\\ -1.086\times10^{+15}\\ +1.086\times10^{+15} $$
Now, reading through Mathematica's documentation, I was expecting the custom made code and MachinePrecision to agree since
The typical arrangement is that all machine‐precision numbers in the Wolfram Language are represented as "double‐precision floating‐point numbers" in the underlying computer system. On most current computers, such numbers contain a total of 64 binary bits, typically yielding 16 decimal digits of mantissa.
And I understand that due to a mantissa of around 16 digits, the first numerical problems would start at around $0.01$, since the biggest number is of order $10^{15}$. This expectation is confirmed by the differences between the numbers obtained within Mathematica with different precisions. But I cannot make sense of differences between the C code and MachinePrecision numbers obtained in Mathematica.
My questions are then:
What numbers should I trust?
Why is MachinePrecision result different than C-code result?
Notes:
a) 'precision' is used as a defining variable whenever I can set Precision or WorkingPrecision. This includes the RandomReal[] functions calls, N[], etc. A numerical generated matrix $M$ has
Precision[M]=precision