Any reliable polygon normal calculation code?

Do you have any reliable 3d polygonal face normal calculation code?

Any language will do, I'll port it to make it work. Even if you find some code in a 3d game engine and post it here I'll be more than grateful.

Edit:

I'm using this but it fails when faces are 90 degrees upright or similar.

        // the normal point
        var x:Number = 0;
        var y:Number = 0;

        var z:Number = 0;

        // if is a triangle with 3 points
        if (points.length == 3) {

            // read vertices of triangle
            var Ax:Number, Bx:Number, Cx:Number;
            var Ay:Number, By:Number, Cy:Number;
            var Az:Number, Bz:Number, Cz:Number;
            Ax = points[0].x;   Bx = points[1].x;   Cx = points[2].x;

            Ay = points[0].y;   By = points[1].y;   Cy = points[2].y;
            Az = points[0].z;   Bz = points[1].z;   Cz = points[2].z;

            // calculate normal of a triangle
            x = (By - Ay) * (Cz - Az) - (Bz - Az) * (Cy - Ay);
            y = (Bz - Az) * (Cx - Ax) - (Bx - Ax) * (Cz - Az);
            z = (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax);

        // if is a polygon with 4+ points
        }else if (points.length > 3){


            // calculate normal of a polygon using all points
            var n:int = points.length;          
            x = 0;
            y = 0;
            z = 0

            // ensure all points above 0
            var minx:Number = 0, miny:Number = 0, minz:Number = 0;
            for (var p:int = 0, pl:int = points.length; p < pl; p++) {

                var po:_Point3D = points[p] = points[p].clone();
                if (po.x < minx) {  minx = po.x;     }
                if (po.y < miny) {  miny = po.y;     }
                if (po.z < minz) {  minz = po.z;     }
            }
            if (minx > 0 || miny > 0 || minz > 0){
                for (p = 0; p < pl; p++) {
                    po = points[p];
                    po.x -= minx;
                    po.y -= miny;

                    po.z -= minz;
                }
            }

            var cur:int = 1, prev:int = 0, next:int = 2;
            for (var i:int = 1; i <= n; i++) {

                // using Newell method
                x += points[cur].y * (points[next].z - points[prev].z);
                y += points[cur].z * (points[next].x - points[prev].x);

                z += points[cur].x * (points[next].y - points[prev].y);
                cur = (cur+1) % n;
                next = (next+1) % n;
                prev = (prev+1) % n;
            }
        }

        // length of the normal
        var length:Number = Math.sqrt(x * x + y * y + z * z);


        // if area is 0
        if (length == 0) {
            return null;

        }else{
            // turn large values into a unit vector
            x = x / length;
            y = y / length;
            z = z / length;
        }

Edit:

I can somehow display 3d arrows for each calculated normal and go so far as to document test cases where it worked and failed. But I was looking for a simpler solution, simply by using code that someone else has created and found to be working. Is OpenGL all you have? Doesn't anyone have any triangle/polygonal normal calculation code at all?

To Jari:

I'm using this code to calculate normals, as built from OpenGL's pseudo code. Is there something I've done wrong? and what would happen if the points are not perfectly inplane? would it fail and return a random normal or just change the normal slightly?

            n = points.length;
            x = y = z = 0;
            for (i = 0; i < n; i++) {
                j = (i + 1) % n;
                x += (points[i].y - points[j].y) * (points[i].z + points[j].z);

                y += (points[i].z - points[j].z) * (points[i].x + points[j].x);
                z += (points[i].x - points[j].x) * (points[i].y + points[j].y);
            }

Answer

I'm successfully using this for quads / triangles.

Loop through each face, and pass in 3 verticies. If you have a quad ABCD pass in ABD. For example for the front facing face on this cube, I would pass in, v2, v3, v0

// cube ///////////////////////////////////////////////////////////////////////
//    v6----- v5
//   /|      /|
//  v1------v0|

//  | |     | |
//  | |v7---|-|v4
//  |/      |/
//  v2------v3

This is in C++ which supports operator overloading, so if in AS3 convert to Vector3d, which I notice you're not using now, might as well. ( I could see why you aren't, sometimes when trying to make found code work i also want to rule out any variables that could make it fail )

//  Modified from http://www.fullonsoftware.co.uk/snippets/content/Math_-_Calculating_Face_Normals.pdf
Vec3f RibbonMesh::calcNormal( const Vec3f &p1, const Vec3f &p2, const Vec3f &p3 )
{
    Vec3f V1= (p2 - p1);

    Vec3f V2 = (p3 - p1);
    Vec3f surfaceNormal;
    surfaceNormal.x = (V1.y*V2.z) - (V1.z-V2.y);
    surfaceNormal.y = - ( (V2.z * V1.x) - (V2.x * V1.z) );
    surfaceNormal.z = (V1.x*V2.y) - (V1.y*V2.x);


    // Dont forget to normalize if needed
    return surfaceNormal;
}

Edit: Found source where I got this from Edit: corrected calculation of z component of surface vector