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](http://mathematica.stackexchange.com/users/15996/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
*)