First some references that are somewhat related. In all of them the large degree of the polynomials is a factor, in addition to the large coefficients. In this question, only the large coefficients are relevant.
On large and small residuals
When measuring the size of a residual, one should consider what scale is appropriate. For machine-precision work, the inherent relative error is usually about 10^-16 (since IEEE double-precision reals are the common floating-point numbers used by Mathematica). So when one finds a root $x = x0$ of $f(x) = 0$,
a residual of about $\big|\,f(x0)\,\big| \approx \big|\,f'(x0)\,\big|\, \times \, 10^{-16}$ would not be uncommon, assuming the root $x0$ was successfully computed and all its bits correct. In a numerical procedure there may also be rounding or other sources of error, that means only most of the bits are correct. The best one could hope for is that $\big|\,f(x0)/f'(x0)\,\big|$ is not more than $10^{-16}$ or so.
The following shows the relative error of the computed roots. There are about 400 that are exactly 0. (the RealExponent of 0. is Log10[$MinMachineNumber] or -307.653). And there seems to be a large group that are off by 1 bit or less. Finally, there is a third group which have a more significant error.
PolyPlot[-400., 400, 1.];
errordata =
Table[poly/Norm@Grad[poly, {x, y}] /. Thread[{x, y} -> xy0], {xy0,
pairs}] // Abs // RealExponent;
MinMax[errordata]
ListPlot[errordata]
(* {-307.653, -12.6089} *)
So most of the solutions returned by Solve are highly accurate.
Count[errordata, x_ /; x > -16]
Count[errordata, x_ /; x <= -16]
(*
1243
1961
*)
On the proper use of SetPrecision
I hope you will forgive me for preaching a little bit.
I have some disdain for the use of SetPrecision[expr, Infinity] (or increased precision) when it is for the purpose showing that the apparent error can be made to disappear by converting the problem to an exact problem, which exact solvers will solve exactly.
The original problem is given in approximate form, and I assume that means the numbers were derived from data and are only approximately known.
(Indeed, one should be careful about assuming the double-precision coefficients are as accurate as their precision. In fact, Mathematica has ways to deal with the known accuracy of data, but I rarely see it used on this site.)
Converting an approximate to an exact problem might make one feel that Mathematica magically found the exact solution.
If so, one has been taken in.
The uncertainty in inputs of a problem should be allowed to propagate through the computed solution and show up in the answers and checks. This is what I tried to
show above.
Increasing the precision of the expression(s) in a problem seems to me to be falsely claiming to know more about the inputs than one does. It seems better to me to understand the results and the underlying numerics, especially in this cases, where, I would claim, the results are correct.
On the other hand, I think it is okay to set the precision to infinity to analyze an approximate problem, as I did in the second linked answer above. Another exception is when the original problem had exact coefficients to begin with (say, from a mathematical problem, not a scientific one). In that case, one should probably not have converted them to machine-precision numbers.
Another way to use SetPrecision[] to analyze this problem is to use SetPrecision[expr, $MachinePrecision]. This converts the machine-precision coefficients to arbitrary-precision numbers that are declared to be accurate to machine precision. It also turns on Mathematica's built-in precision tracking, so that one can estimate the error by examining the computed results.
Below we take the MachinePrecision results and check the residuals after converting them to $MachinePrecision. (Note MachinePrecision and $MachinePrecision are not the same.)
residualsMP =
Table[SetPrecision[poly, $MachinePrecision] /.
Thread[{x, y} -> xy0], {xy0, SetPrecision[pairs, $MachinePrecision]}] // Abs;
residualsMP // Short[#, 5] &
AllTrue[residualsMP, PossibleZeroQ]
(* True *)
They all turn out to be "zero." In fact, the form "0. x 10^16" means that the rounding error is around 10^16 and the computed value is closer to zero than that. We can examine these in closer detail. Let's look at the first and third residuals, which have different results displayed above.
residualsMP[[{1, 3}]] // InputForm
(*
{0``-16.998893852353977,
0``-5.4477916593052464}
*)
The computed (approximate) errors are 10^16.998893852353977 and 10^5.4477916593052464 (the negative sign is just how Accuracy works).
These should be bounds on the computed solutions. As we can see below, the computed machine-precision residuals are less than that:
poly /. Thread[{x, y} -> pairs[[1]]] // Abs // Log10
poly /. Thread[{x, y} -> pairs[[3]]] // Abs // Log10
(*
16.1307
3.86662
*)
