I'm getting WA #12 in this problem, what I do is take to points that are consecutive in the convex hull, then order all the other (N-2) points according to the measure of the angle we get of joining this point with the other two I chose at the beginning.

#include <iostream>
#include <cmath>
#include <algorithm>

using namespace std;

struct point{
    long long x,y;
    long double ang;
    
    point(){}
    
    bool operator < (point X) const{
        return ang < X.ang;
    }
}P[5000];

bool ccw(point &A, point &B, point &C){
    return (A.x * B.y + B.x * C.y + C.x * A.y - A.y * B.x - B.y * C.x - C.y * A.x) >= 0;
}

long double dist(point &A, point &B){
    return sqrt((long double)(A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
}

int main(){
    int N;
    
    cin >> N;
    
    for(int i = 0;i < N;++i)
        cin >> P[i].x >> P[i].y;
    
    for(int i = 1;i < N;++i)
        if(P[i].x < P[0].x)
            swap(P[0],P[i]);
    
    for(int i = 2;i < N;++i)
        if(ccw(P[0],P[1],P[i]))
            swap(P[1],P[i]);
    
    long double a,b,c = dist(P[0],P[1]);
    
    for(int i = 2;i < N;++i){
        a = dist(P[i],P[0]);
        b = dist(P[i],P[1]);
        
        P[i].ang = (a*a + b*b - c*c) / (2 * a * b);
    }
    
    sort(P + 2,P + N);
    
    cout << P[0].x << " " << P[0].y << endl;
    cout << P[1].x << " " << P[1].y << endl;
    cout << P[(N + 1) / 2].x << " " << P[(N + 1) / 2].y << endl;
    
    return 0;
}