This is my source : I think everything is correct!!!

#include <fstream.h>
#include <stdlib.h>
#include <math.h>
#include <iomanip.h>

  struct TPoint
    {double x, y;};

  struct TAngle
    {
      double ang;
      int ind;
    };

  TPoint A[5100];
  TAngle B[5100];
  int n;
  double min;
  int ind;
  TPoint tmp;

  void ReadData ()
    {
      cin >> n;
      for (int i = 0; i <= n - 1; i++)
        {cin >> A[i].x >> A[i].y;}
    }

  double dist (TPoint const &a, TPoint const &b)
    {return sqrt((a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y));}

  int comp (const void *e1, const void *e2)
    {
      double z = ((TAngle*)e1)->ang - ((TAngle*)e2)->ang;
      if (z < 0)
        return -1;
      if (z > 0)
        return 1;
      if (z == 0)
        while (1) // Here if there is four points on one circle
program gets TL...
          {}
    }

  void Solve ()
    {
      min = A[0].y;
      ind = 0;
      for (int i = 1; i <= n - 1; i++)
        if (A[i].y < min)
          {
            min = A[i].y;
            ind = i;
          }
      tmp = A[0];
      A[0] = A[ind];
      A[ind] = tmp;
      min = atan2 (A[1].y - A[0].y, A[1].x - A[0].x)*1000;
      ind = 1;
      for (int i = 2; i <= n - 1; i++)
        if (atan2 (A[i].y - A[0].y, A[i].x - A[0].x)*1000 < min)
          {
            min = atan2 (A[i].y - A[0].y, A[i].x - A[0].x)*1000;
            ind = i;
          }
      tmp = A[1];
      A[1] = A[ind];
      A[ind] = tmp;
      double a, b, c;
      for (int i = 2; i <= n - 1; i++)
        {
          B[i].ind = i;
          a = dist (A[0], A[i]);
          b = dist (A[1], A[i]);
          c = dist (A[0], A[1]);
          B[i].ang = acos ((a*a + b*b - c*c)/(2*a*b));
          B[i].ang *= 1000;
        }
      qsort (&B[2], n - 2, sizeof (TAngle), comp);
    }

  void WriteData ()
    {
      cout << setiosflags (ios::fixed) << setprecision (0) << A[0].x
<< " " << A[0].y << endl;
      cout << setiosflags (ios::fixed) << setprecision (0) << A[1].x
<< " " << A[1].y << endl;
      cout << setiosflags (ios::fixed) << setprecision (0) << A[B[(1
+ n)/2].ind].x << " " << A[B[(1 + n)/2].ind].y << endl;
    }

  int main ()
    {
      ReadData ();
      Solve ();
      WriteData ();
      return 0;
    }