Dealing with interpolated normals within a raytracer
Disclaimer: It's been a long time since I looked at this sort of thing but here goes...
Disclaimer2: On re-reading your question(s) I realised I might have misunderstood what you were asking. I'll leave this here just in case you were looking for this sort of reply.

Are you rendering from the original Bezier patch definition (e.g. https://www.cs.utah.edu/gdc/projects/alpha1/help/man/html/model_repo/model_teapot/teapot.pts)  or from a (post tessellation) triangle representation?

If the former, then you should, in the vast majority of cases, be able to 'trivially' generate the normals directly from the Bezier representation, by computing surface partial derivatives at each ray intersection point, to give tangent vectors, and computing the cross product. The cases you do have to watch out for are at very centres of the bottom and top of the teapot where the patches, though geometrically well-defined, have a zero partial derivative because the rectangular patch has been squeezed into a triangle, e.g. control points for U,V = (0,0), (0,1/3), (0,2/3), (0,1) are all identical, so dPos/dv is 0.  All you need do is compute dPos/du at (0,0) and (0,1) to get valid tangents.

If you are rendering from triangles, I assume you already have normals at the triangle vertices. If you use an appropriate ray-triangle intersection test you can get the barycentric coordinates of the intersection point. Use these to blend the vertex normals (and after renormalising) you'll have a 'reasonable' surface normal (though it will only be C0 continuous).