OK, here s a stab at a general algorithm to carry out this task, independent of what language and graphics library you are using.

1. I am going to assume you want to calculate vertex normals i.e. you will end up with 1 normal per vertex
2. As well as your 65000 vertex positions, I am also assuming you have some list of indices which define the triangles in the mesh e.g. triangle 1 in the mesh is made of vertices 3, 7 & 24. Whether you have a triangle list or strip or some other description of the meshes faces, they all basically amount to the same thing: they describe some sort of method for determining which vertices are in a particular face. The difference between such methods is often relating to efficiency, but it does not really change the basic idea behind this algorithm.
3. The first thing you need to do is calculate the normal for each triangle in the mesh. For a triangle with vertex position vectors a, b & c, you calculate the edge vectors for 2 of the edges in the triangle e.g. edge1 = b - a, edge2 = c - a. You then take the cross product of these 2 vectors to get the normal vector to this triangle. The result of the cross product is a vector that will need to be normalised.
4. Once all triangle normals have been calculated, the normal for a vertex is calculated by averaging the normals for each triangle that the vertex belongs to. You can just do a straight unweighted average of the normals i.e. if a vertex is part of 3 faces with normals n1, n2 & n3, the vertex normal is simply (n1 + n2 + n3)/3. You can also do a weighted average, where each of the triangle normals in the sum are weighted by some importance factor like the area of that triangle. There s no correct way here and you can play around with different things. Either way, once the average has been calculated, it also needs to be normalised. This however can be done in 1 step since, for 3 vectors n1, n2 & n3, normalised(n1 + n2 + n3) == normalised((n1 + n2 + n3)/3).

Now, I must stress that this is a rough description of what needs to be accomplished in order to do this there are certain corners that can be cut for efficiency. What I m getting at is that you don t need to calculate all triangle normals first and then calculate all vertex normals - the 2 steps can be mixed in to make things more efficient.

Some pseudo code might look this

Vector verts[65000];         // 65000 vertex positions
Triangle faces[87000];       // as an example, 87000 triangles
Vector normals[65000];       // 1 normal per vertex - initialised to (0, 0, 0)

// loop over all faces to calculate vertex normals
for (int i=0 ; i<87000 ; i++)
{
    Vector v1 = verts[faces[i].vertex1];
    Vector v2 = verts[faces[i].vertex2];
    Vector v3 = verts[faces[i].vertex3];

    Vector edge1 = v2 - v1;
    Vector edge2 = v3 - v1;
    Vector normal = edge1.CrossProduct(edge2);  // or edge2.CrossProduct(edge1)
    normal.Normalise();

    normals[faces[i].vertex1] += normal;
    normals[faces[i].vertex2] += normal;
    normals[faces[i].vertex3] += normal;
}

// vertex normals need to be normalised
for (int i=0 ; i<65000 ; i++)
{
    normals.Normalise();
}

Regarding some of the other comments about winding order - if you get this wrong, the normals will point inwards, rather than outwards. This can be fixed by simply changing the order of the cross product, as noted above.

since you have the indices, i am assuming its either a triangle list / strip or fan. Read each triangle. Compute the normal by taking the cross product of 2 of the vectors of the triangle. You have 1 problem though. If you dont know the winding order of the triangles, then you might get the opposite value. Which software created the mesh ? Can you inspect within the data file or software what the winding order was ? is it left or right handed ?

What you need is simply the cross-product of the vectors that make up two sides of the triangle, normalized to a unit vector.

As Andrew Keith said in his comment, you d better know that the triangles are defined either clockwise or counterclockwise when you look at the triangle from the "outside." If you can t guarantee consistency, you have a mess on your hands. But probably (or at least hopefully) the code that created the object is sane.

Steg s answer points into the right direction. However, if you need high-quality normals, have a look at the paper Discrete Differential-Geometry Operators for Triangulated 2-Manifolds. The cotangent formula (8) gives you good results even for irregular meshes where estimates such as the triangle area break down.