One way to do this would be via a vertex weight proximity technique.
Start by scaling plane up a billion-fold or whatever (about the median of its vertices would be fine), so that the face actually represents the plane, rather than the bounded plane (face.) Then, give the cube a new vertex group and assign all verts to it. Give the cube a vertex weight proximity modifier, targeting plane, referencing that vertex group, on geometry/face mode, with a highest of 0.0 and a lowest of 1E-31 (or whatever). Apply the modifier, then run a weight paint/clean operation if needed. Only verts that were within 1E-31 units of the plane remain assigned to that vertex group.
Clearly, with the x1 billion scale and the 1E-31 proximity modifier, this is valid only to a certain precision, but you can get any particular precision you'd like (and we're talking floating point values, so there's not really such thing as perfect precision anyways.)
I'm not sure that I understand your "even better" question, but if A, B and C are points on cube that you want to select on the basis of proximity to the triangular prism, I think you can see that you can use the same technique (obviously, don't scale the prism up a billion-fold.)
Not Python, but you indicated an openness to any technique, and of course all of the operations I mentioned should be perfectly scriptable.
Just to make sure you know, that the fact that a point exists on a plane does not mean that any faces that it is part of are co-planar with that plane. A might be on the plane while B, C and D are not.
(And if you're curious, how I'd handle the linked problem is by running a shrinkwrap->outside followed a vertex weight prox. I don't think it's the kind of answer the original asker was looking for, but it may be the kind of answer you're looking for.)