Each vertex belongs to one or more faces (usually triangles, sometimes quads -- I'll use triangles in this answer).
A triangle that is not attached to any other triangles cannot be 'smoothed'. It is flat. Only when a face has neighbours can you reason about smoothing them together.
For a vertex where multiple faces meet, calculate the normals for each of these faces. The cross product of two vectors returns a perpendicular (normal) vector, which is what we want.
A --- B
\ /
C
v1 = B - A
v2 = C - A
normal = v1 cross v2
Be careful to calculate these vectors consistently across all faces, otherwise your normal may be in a negative direction to that you require.
So at a vertex where multiple faces meet, sum the normals of the faces, normalise the resulting vector, and apply it to the vertex.
Sometimes you have a mesh where some parts of it are to be smoothed, and others not. An easy to picture example is a cylinder made of triangles. The round surface of the cylinder would smooth well, but if you consider triangles from the flat ends at the vertices around the sharp ridge, it will look strange. To avoid this, you can introduce a rule that ignore normals from faces which deviate too far from the normal of the face you're calculating for.
EDIT there's a really good video showing technique for calculating Gourad shading, though it doesn't discuss an actual algorithm.
You might like to take a look at the source of of Three.js. Specifically, the computeVertexNormals
function. It does not support maintaining sharp edges. The efficiency of your algorithm depends to a large extent upon the way in which you are modelling your primitives.