well like I said, GLOBAL, if the user doesn't want to sit and wait for eternity for this to calculate they should have the option of using a less accurate method, such as just providing an angle to use on all polygons.
norm -1, I would bet rounding error is responsible for this, in stead of using vector==vector, you should use dot(vector,vector)>0.99f or something similar.
I have a geometry class in the works that is built around making a lot of these comparisons faster, it has a single list of points which is referenced by a list of verts which are used by a list of polygons, every point keeps a reference to every vertex which uses it, and every vertex keeps a reference to every polygon that uses it, it makes looking up stuff like neighbors extremely fast. it'l be (nearly) an N*logN operation rather than (worse than) N^2, because the largest chunck of data (points) will work on a ordered index, and you can use binary search on that. (keeping everything in order is the major coding problem of this method, especially given the interdependencies, but it's the only way to avoid an unacceptably long processing time). I have it set up so it can keep everything up to date all the time (this mode is mostly for geometry editing) or allow you to make a few changes and then fix it later.