Return to Answer

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.

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.

[Update: D'oh -- sometimes I get so focused on something of interest in the OP's approach I can't see the forest for the trees. I should have thought, we can use Solve, as shown below, because it's a polynomial...]

##Two explanations##

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];

With the other roots, it is more or less the same, that (1) the computed number for x is the best or nearly best possible approximation of the root and (2) the value of e is due to a round-off error and/or a large value of D[e, x]. 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 cases are the complex conjugate pair sol[[{57, 58}]] which are 12 ϵ away from the best solution (as found by FindRoot). Below we see the value of e is less at x - 12 eps;

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;

[Update: D'oh -- sometimes I get so focused on something of interest in the OP's approach I can't see the forest for the trees. I should have thought, we can use Solve, as shown below, because it's a polynomial...]

##Two Three explanations##

There are two three 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];

###Update: Using Solve###

If we use Solve on the exact equation, we can check these solutions against the ones returned by NSolve. The equation has to have exact coefficients or Solve uses numerical techniques.

solexact = Solve[SetPrecision[e, Infinity] == 0, x];

(x /. sol) == Sort@N[x /. solexact]
Sort@N[x /. solexact] - (x /. sol) // Abs // Max
(*
  True
  3.55271*10^-15
*)

Here we see Solve did in fact return solutions (of course!).

ParallelMap[FullSimplify, 
  SetPrecision[e, Infinity] /. solexact] // AbsoluteTiming
(*
  {77.192719, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
     0, 0, 0}}
*)

Perhaps the answer should stop here. Originally I was interested in explaining the numerics because this issue has come up before. I thought having an explanation of why NSolve appears to have done a bad job but in fact did a good job might be helpful.

With the other roots, it is more or less the same, that (1) the computed number for x is the best or nearly best possible approximation of the root and (2) the value of e is due to a round-off error and/or a large value of D[e, x]. 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 cases are the complex conjugate pair sol[[{57, 58}]] which are 12 ϵ away from the best exact solution (as found by FindRoot Solve). Below we see the value of e is less at x - 12 eps;

dxeps = MapThread[
  Function[{x, xexact}, 
   Through[{Re, Im}[(xexact) - (x)]] /  
      {2^Floor@Log2@Abs[$MachineEpsilon Re[x0]], 
       2^Floor@Log2@Abs[$MachineEpsilon Im[x0]]} /. x0 -> x
   ], 
  {x /. sol, Sort@N[x /. solexact]}
  ]

Update: Perhaps one should just use Solve as above. I'll leave the explanation below because I find it interesting, even if it's rather useless in this case.