I would pick any point on the shape* as starting value of extremes and iterate over all other points on shape testing whether they lie on left/right side of line from P to left/right extreme respectively, reassigning value if new extreme is found:
void extreme_points(Point& left, Point& right, Point& p, std::vector<Point> polygon)
{
left = right = polygon[0];
for(auto v : polygon)
{
if(isLeft(p, left, v) > 0.0)
{
if(abs(isLeft(p, right, v)) > Epsilon)
left = p;
else { /*handle special case where P is parallel,*/}
}
if(isLeft(p, right, v) < 0.0)
{
if(abs(isLeft(p, right, v)) > Epsilon)
right = p;
else { /*handle special case where P is parallel,*/}
}
}
}
where isLeft is well-known "sign test". Example of implementation:
float isLeft( Point P0, Point P1, Point P2 )
{
return ( (P1.x - P0.x) * (P2.y - P0.y)
- (P2.x - P0.x) * (P1.y - P0.y) );
}
You may also want to handle cases for almost zero values(all three points are on the same line) picking, for example, the closer one. You can also find a library with more accurate isLeft implementation.
*this should work for any convex and star-shaped polygons.