Why do 3d modeling packages use a single normal per polygon and how is this viable for smooth shading
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.