--- edit ---
Okay, I missed this subtlety of "infinite" solutions. They are indeed infinite, in an important sense (and this gives the numeric instability). In a certain way the system is not "generic". Specifically, coefficients satisfy one or more (unknown to me) relevant algebraic conditions that put the system on what is called the "discriminant variety" for the family of systems having the same Newton polytope. Remark: Since the system is generic for your purposes, that means your entire family of systems lives on this variety.
In such a situation the number of solutions can drop from the generic (in this case) 8. But that does not quite happen for the example above, because we have coefficients given only approximately, via machine doubles. So we actually have a system NEAR but not ON the discriminant variety.
So what happens to the 8 solutions when we are on this variety? In this case 6 have, loosely speaking, wandered off to infinity. Continuity of solutions as a function of coefficients implies that systems near the variety will have 6 "large" solutions. That's what NSolve is showing. When I round and solve to high precision I get a perfectly good solution set, albeit residuals of the "bad" solutions will be orders of magnitude larger than those of the good ones.
Homotopy continuation methods will not show behavior this for either of two reasons.
(1) They detect "diverging" solutions via heuristics in the tracking phase. In this case 6 would be dropped once they got too large in the tracking.
(2) In some cases one already has a correct result for a related system from the same part of the discriminant variety. In this case a common tactic is to use what is called the "Cheater's Homotopy", using these old solutions to go to solutions for the new system. Here we know in advance that we have the right number of starting points and the tracking will work just fine (barring other numeric issues).
NSolve has a system option that can be used, in a heuristic manner, to elicit better behavior. Or worse, depending on mood. The idea is implementation dependent: we compute a numerical Groebner basis, and a nondefault setting will tell that code to set to zero any coefficients that are below some threshold smaller than the average size in the rest of the polynomial. I show how to use this below, with 10^(-6) for that tolerance. Remark: if I go to 10^(-4) it will hang. That's one way the heuristic can mess up. Another is that it can give garbage results. It is very much a situation where one needs to make sure the results (when they appear) seem plausible. Even finding a viable tolerance is a trial-and-error process, informed perhaps by the nature of the problem. Here I know we have approximated coefficients (only) to several digits, so I expect that 10^(-8) or so will work for my purpose. And it did, as did 10^(-6).
SetSystemOptions["NSolveOptions" -> {"Tolerance" -> 10^(-6)}];
sol = NSolve[e == 0, var]
(* {{c[1] -> -0.840799, c[2] -> 0.0280179, c[3] -> 0.877839,
s[1] -> 0.541347, s[2] -> -0.999607, s[3] -> -0.478956,
y1[1] -> 17.8457, y1[2] -> 0.56958, y1[3] -> -0.594672,
y2[1] -> -0.594672, y2[2] -> 17.8457,
y2[3] -> 0.56958}, {c[1] -> 0.840799, c[2] -> -0.0280179,
c[3] -> -0.877839, s[1] -> -0.541347, s[2] -> 0.999607,
s[3] -> 0.478956, y1[1] -> -17.8457, y1[2] -> -0.56958,
y1[3] -> 0.594672, y2[1] -> 0.594672, y2[2] -> -17.8457,
y2[3] -> -0.56958}} *)
Here is the example where I first ran into this near-discriminant-variety computational mishap four years ago.
http://homepages.math.uic.edu/~jan/
It took some head-scratching between myself and the web page author to figure out what had gone wrong.
--- end edit ---
