If the constraint is that all of the output points be on the surface, you want a consistent method of addressing the surface itself rather than worrying about the 3d > surface conversion for your points.
The hacktastic way to do that would be to create a UV map for your 3d object, and then scatter points randomly in 2 dimensions (throwing away points which happened not to land inside a valid UV shell). Once your UV shells are filled up as much as you'd like, you can convert your UV points to barycentric coordinates to convert those 2-d points back to 3-d points: effectively you say "i am 30% vertex A, 30 % vertex B, and 40% vertex C, so my position is (.3A + .3B + .4C)
Besides simplicity, another advantage of using is UV map is that it would allow you to customize the density and relative importance of different parts of the mesh: a larger UV face will get a lot of scattered points, and a smaller one fewer -- even if that doesn't match the physical size or the faces.
Going to 2D will introduce some artifacts because you probably will not be able to come up with a UV map that is both stretch-free and seam-free, so you'll get variations in the density of your scatter because of that. However for many applications this will be fine, since the algorithm is really simple and the results easy to hand tune.
I have not used this one but this looks like it's based on this general approach: http://www.shanemarks.co.za/uncategorized/uv-scatter-script/
If you need a more mathematically rigorous method, you'd need a fancier method of mesh parameterization : a way to turn your 3-d collection of triangles into a consistent space. There is a lot of interesting work in that field but it would be hard to pick a particular path without knowing the application.