Basically, there is nothing wrong with the result that Solve/NSolve produces, or at least the errors are slight. The number 10^166 may not seem that close to zero, but that depends on whether it is compared with 1 or 10^175. Since NSolve solves for x, the best you can hope for is that all the bits in the machine number for x are correct. In practice, the last few bits may be wrong due to round-off error in computing function values.
Let me review a few facts about machine floating-point numbers that I wish to refer to in my explanation. A (binary) machine real x is a rational number of the form x0 = mantissa * 2^exponent, where mantissa is represented with a fixed number of bits. Changing the lowest-order bit corresponds to changing x by the smallest amount ϵ possible. The real number x can be thought of as representing a real number that lies between x0 ± ϵ/2, with a maximum uncertainty of ϵ/2; the point-estimate x0 is the number used in computation. In a computational process there can be round-off error that gets propagated resulting in a number x that has a uncertainty greater than ϵ/2. In Mathematica, the ϵ for the number x = 1 is given by $MachineEpsilon (on my machine, $MachineEpsilon == 2.220446049250313*^-16 == 2^(-52)). For a number x, one can estimate the corresponding ϵ with $MachineEpsilon * x. One can get the exact number x0 with x0 = SetPrecision[x, Infinity].
There are two ways to explain why e /. x -> sol[[1]] can result in a number larger than 10^166 when sol[[1]] gives the best possible approximation of the root (at machine precision). For the sake of reproducibility, here is the code I used from the OP, with a slight change (Rationalize did nothing on the coefficients generated):
SeedRandom[1];
(*r := Rationalize[RandomReal[NormalDistribution[0,1]]];*)
r := RandomReal[NormalDistribution[0, 1]];
e = randompol[100, 1];
sol = NSolve[e == 0, x];
First, if y = f[x], then for a given error dx in x, the corresponding error dy in y is approximated by
dy = f'[x] dx
The size of dy evidently depends on the size of f'[x]. If dx is on the order of 0.5 $MachineEpsilon * x and the derivative D[e, x] /. sol[[1]] of the OP's equation is on the order of 10^180, then dy will be on the order of
Abs[D[e, x] (0.5 $MachineEpsilon*x)] /. sol[[1]]
(*
1.77676*10^166
*)
The magnitude of e, which is the error from 0, is about twice that,
e /. sol[[1]]
(*
-3.82848*10^166
*)
suggesting either round-off error in the calculation of e or error in the result of NSolve. If the error is with NSolve, the error is not much since the uncertainty in x appears to be a little more than $MachineEpsilon * x. It appears that the root in sol[[1]] might be off in its last binary bit from what is optimal. In any case, it's a good-looking result.
But not so fast. A second way to analyze the situation is this. Mathematica can track precision, if so-called arbitrary-precision numbers are used. We can turn precision-tracking on, by setting the precision to $MachinePrecision.
SetPrecision[e, $MachinePrecision] /. SetPrecision[sol[[1]], $MachinePrecision]
Accuracy[%]
(*
0.*10^168
-168.799
*)
Mathematica is indicating that the result is 0 within an error less than about ±10^168. Our original value -3.8*^166 for e /. sol[[1]] is well within that error. What is happening is a tremendous loss of precision in computing e /. sol[[1]]. If we set the precision high enough, we can calculate the value of e at the (exact rational) number computed for sol[[1]].
SetPrecision[e, 50] /. SetPrecision[sol[[1]], 50];
Precision[%] (* check precision *)
N[%%] (* print as a machine Real *)
(*
31.3568
-1.28981*10^166
*)
By controlling the round-off error in calculating e, we see that the value for e is within the estimated maximum error 1.77676 * 10^166. This suggests that the computed value for x is the best possible. We can check that with the number within ±ϵ. They give values of e further from zero:
With[{eps = SetPrecision[$MachineEpsilon, 50],
sol = SetPrecision[sol[[1]], 50]},
{SetPrecision[e, 50] /. x -> (x - eps x /. sol),
SetPrecision[e, 50] /. x -> (x + eps x /. sol)} // N
]
(*
{-4.84333*10^166, 2.2637*10^166}
*)
Thus we see that NSolve actually computed the best approximation to the real root for the case sol[[1]]. It is more or less the same with the other roots. About half give the best possible solution. In most other cases, the error is in the last one or two bits, but that might not be surprising given the sort of numerical difficulty e presents. The worst case are the conjugate pair sol[[{57, 58}]] which are 12 ϵ away from the best solution (as found by FindRoot).
e /. sol[[57]] // Abs
With[{eps = 2^Floor@Log2@Abs[$MachineEpsilon Re[x]] /. sol[[57]]},
e /. x -> (x - 12 eps /. sol[[57]])
] // Abs
(*
3.16743*10^46
6.19979*10^45
*)
The code finds how far each solution from NSolve is from each root.
dxeps = Function[{sol},
Through[{Re, Im}[(x /. #) - (x /. sol)]] /
{2^Floor@Log2@Abs[$MachineEpsilon Re[x]],
2^Floor@Log2@Abs[$MachineEpsilon Im[x]]} /. sol &@
FindRoot[SetPrecision[e, Infinity], {x, x /. sol},
WorkingPrecision -> 200]
] /@ sol;
48/100 are best possible; the worst case are solutions 57/58.
Count[dxeps, {0., 0.}]
(*
48
*)
Position[#, Max@#] &@Abs@dxeps
(*
{{57, 1}, {58, 1}}
*)
Distribution of errors
Labeled[
Histogram[Flatten@Abs@dxeps],
Style["Errors in real and imaginary parts of roots (epsilons)", "Label"]
]
As mentioned by user15996, to get solutions that produce values for e that are zero to some number of digits of accuracy, greater working precision is needed. First, we will need to set the precision of the coefficients of the polynomial e to be exact (or of sufficient precision), which I'll call eP. As an example let's consider sol[[57]]. We see below that there is a deficit of almost 49 digits, that is, Mathematica is telling us that the value is zero with an error less than 10^49. So if we use a WorkingPrecision of 49 digits plus $MachinePrecision, we should end up with an accuracy around 0 digits. It's a little tricky to get exactly 0 because the real and imaginary parts have different accuracies. We can increase the WorkingPrecision to get a positive number of significant digtis.
eP = SetPrecision[e, Infinity];
eP /. SetPrecision[sol[[57]], $MachinePrecision]
accuracydeficit = Max[-Accuracy /@ {Re[%], Im[%]}]
(*
0``-48.9342823110845 + 0``-48.95998064695299 I
48.96
*)
extraprecision = 0;
extraprecision + accuracydeficit + $MachinePrecision
solWP = NSolve[eP == 0, x,
WorkingPrecision -> extraprecision + accuracydeficit + $MachinePrecision];
eP /. solWP[[57]]
(*
64.9146
0``0.025698335868500944 + 0``0 I
*)
extraprecision = 2;
solWP = NSolve[eP == 0, x,
WorkingPrecision -> extraprecision + accuracydeficit + $MachinePrecision];
eP /. solWP[[57]]
(*
0``2.025698335868505 + 0``2. I
*)
The value of the polynomial with MachinePrecision coefficients is still way off, for the same reasons as before, round-off error with insufficient precision.
e /. solWP[[57]]
(*
-3.37129*10^45 - 2.92341*10^45 I
*)
Therefore to get solutions that will evaluate eP to $MachinePrecision accurately, we should pick the maximum accuracy deficit:
accuracydeficit =
Max[-Map[Accuracy,
Through[{Re, Im}[eP /. SetPrecision[sol, $MachinePrecision]]], {-1}]]
extraprecision = $MachinePrecision;
extraprecision + accuracydeficit + $MachinePrecision
solWP = NSolve[eP == 0, x,
WorkingPrecision -> extraprecision + accuracydeficit + $MachinePrecision];
(*
168.795
200.704
*)
Check:
Through[{Re, Im}[eP /. solWP]] // Max
(*
0.*10^-16
*)