Why do 3d modeling packages use a single normal per polygon and how is this viable for smooth shading If we have a quad that is non-planar. Why do most 3d modeling packages choose to use a single normal for the quad (Specifically for Blender3d by using cross product of two edge's vectors that do not share a vertex, here's the math). Instead of viewing it as two triangles, and using a separate normal for each triangle?
Normally we calculate a vertex normal as the normalized average of all polygons the vertex is part of. If we have a quad consisting of two triangles, then the two vertices that are part of both triangles, get the average of the two triangle normals. This results in smooth shading across the non-planar quad. Instead of the uniform shading we would get if we were to use a single normal for the whole quad.
This also extends to polygons in general. In 3d modeling sofware a single normal is calculated for the whole polygon (for Blender3d with Newell's method), instead of individual normals for each triangle it consists of.
To build on top of this, if we have a mesh with flat shading, each polygon gets a single normal as described above. Flat shading means that the shared vertices are split. And each vertex gets the normal of the polygons it's part of. However what happens if we want to apply a smooth shading on the mesh? Are the vertex normals calculated as average of all polygons it is part of? This doesn't seem right. Do we instead triangulate the mesh, and calculate the vertices as the average of all triangles it's part of?
I think the problem comes from this starting assumption:

If we have a quad that is non-planar.

A non-planar quad does not have a normal. It's not a flat surface, so you can't talk about what its normal would be. You can't talk about the normal of a sphere, a cone, or any other non-flat surface. You can talk about what the normal might be at any particular point on that surface, but there is no normal for the surface in its entirety.
You are asking a question for which there is no correct answer. But there are many incorrect answers, and they're all equally wrong. So there's no reason to pick one wrong answer over another.
If you're talking about computing vertex normals, a non-planar quad is still not a valid surface. The two triangulations for a non-planar quad represent very different surfaces, which would have very different vertex normals. Neither answer can claim legitimacy over the other, so both answers are equally wrong.
So why choose that particular method? Because it's no more wrong than any other method you might pick.
It does not make sense to talk about the vertex or face normals for non-planar polygons. So if you want an answer that makes sense, fix that problem first.