Оптимизация кода - C (СИ)

/////////////////////////////////////////////////////////////////////////////////////
//Нужно оптимизировать код, т.к. при сдаче работы есть ограничение по времени - 1с, 
//а к сожалению работа занимает дольше 1 секунды.
//
//Собственно, вот задача
//На плоскости задан набор точек. Определите площадь наибольшего из треугольников, 
//которые можно построить, используя тройки заданных точек в качестве вершин.
//На стандартном потоке ввода задаётся сначала натуральное число N (1 ≤ N ≤ 650), 
//а затем — координаты N точек в виде Xi, Yi. Координаты каждой точки задаются 
//на отдельной строке и разделены пробелом. Все координаты — целые числа, 
//не превосходящие по модулю 100. Точки в наборе могут повторяться.
//На стандартный поток необходимо напечатать единственное число — площадь наибольшего 
//из треугольников, образованных тройками из данного набора точек. 
//Выводите ответ с шестью знаками после десятичной точки.
/////////////////////////////////////////////////////////////////////////////////////
#include <algorithm>
#include <complex>
#include <cstdlib>
#include <ctime>
#include <iostream>
#include <iterator>
#include <limits>
#include <set>
#include <vector>
/////////////////////////////////////////////////////////////////////////////////////
typedef double                 T_coord;
typedef std::complex<T_coord>  T_point;
typedef std::vector<T_point>   T_points;
/////////////////////////////////////////////////////////////////////////////////////
template<class T>
bool  equal_to_for_real(T a, T b) 
{
    const T  coef = 10;
    return abs(a - b) < std::numeric_limits<T>::epsilon() * coef;
}
/////////////////////////////////////////////////////////////////////////////////////
template<class T>
bool  greater_for_real(T a, T b) 
{
    return a > b
           && !equal_to_for_real(a, b);
}
/////////////////////////////////////////////////////////////////////////////////////
template<class T>
bool  less_for_real(T a, T b) 
{
    return a < b
           && !equal_to_for_real(a, b);
}
/////////////////////////////////////////////////////////////////////////////////////
template<class T>
bool  greater_equal_for_real(T a, T b) 
{
    return !less_for_real(a, b);
}
/////////////////////////////////////////////////////////////////////////////////////
template<class T>
bool  less_equal_for_real(T a, T b) 
{
    return !greater_for_real(a, b);
}
/////////////////////////////////////////////////////////////////////////////////////
struct T_point_generator
{
    int mod_;
    //-------------------------------------------------------------------------------
    T_point_generator(int mod) : mod_(mod)
    {}
    //-------------------------------------------------------------------------------
    T_point  operator() ()
    {
        T_coord  X = rand() % mod_;
        T_coord  Y = rand() % mod_;
        return  T_point(X, Y);    
    }
};
/////////////////////////////////////////////////////////////////////////////////////
void  get_points(T_points&  points)
{
    const int MOD = 100;
    std::generate(points.begin(), points.end(), T_point_generator(MOD));
}
/////////////////////////////////////////////////////////////////////////////////////
struct  T_bottom_left_compare
{
    bool operator() (T_point A, T_point B)
    {
        return equal_to_for_real(A.imag(), B.imag()) 
            ? less_for_real(A.real(), B.real()) : less_for_real(A.imag(), B.imag());
    }
};
/////////////////////////////////////////////////////////////////////////////////////
T_point  vect(T_point  A, T_point  B)
{
    return B - A;
}
/////////////////////////////////////////////////////////////////////////////////////
struct  T_polar_corner_compare
{   
    bool  operator() (T_point A, T_point B)
    {
        return equal_to_for_real(arg(A), arg(B)) 
                   ? less_for_real(abs(A), abs(B)) 
                   : less_for_real(arg(A), arg(B));
    }
};
/////////////////////////////////////////////////////////////////////////////////////
T_coord  det(T_point A, T_point B)
{
    return  A.real() * B.imag() - A.imag() * B.real();
}
/////////////////////////////////////////////////////////////////////////////////////
T_points  get_convex_hull_points(const T_points&  points)
{
    //Ищем выпуклую оболочку методом Грэхема (время работы O(n*lg(n))).
    
    //Находим левую из самых нижних точек заданной совокупности.
    T_point  bottom_point 
        = *std::min_element(points.begin(), points.end(), T_bottom_left_compare());    
    
    //Вычитаем из всех точек вектора точку bottom_point (тогда стартовая точка
    //переместится в начало координат) и записываем эти точки в множество, 
    //отсортированное по полярному углу.
    typedef std::set<T_point, T_polar_corner_compare>  T_convex_hull_points;
    T_convex_hull_points  convex_hull_points;
    
    std::transform(points.begin(), points.end(), 
                   std::inserter(convex_hull_points, convex_hull_points.begin()), 
                   std::bind2nd(std::minus<T_point>(), bottom_point));
 
    bottom_point = T_point();
 
    //Помещаем в стек стартовую точку и первую точку множества, удаляя их из множества.    
    T_points  stack;
 
    T_point  cur_point = bottom_point;
    stack.push_back           (cur_point);
    convex_hull_points.erase  (cur_point);
 
    cur_point = *convex_hull_points.begin();
    stack.push_back           (cur_point);
    convex_hull_points.erase  (cur_point);
 
    //В цикле проверяем, какой поворот образуют точки предпоследняя в стеке,
    //последняя в стеке, и текущая. Если это не левый поворот, то выталкиваем
    //точку из стека (если она там не одна), и вставляем туда текущую. 
    //В конце в стеке останутся точки оболочки.
    for(T_convex_hull_points::const_iterator  point_it = convex_hull_points.begin();
        point_it != convex_hull_points.end(); )
    {
        int last_ind       = stack.size()  - 1;
        int prev_last_ind  = last_ind      - 1;
        T_point  v1 = vect(stack[prev_last_ind], stack[last_ind]);    
        T_point  v2 = vect(stack[last_ind], *point_it); 
 
        if(det(v1, v2) <= 0)
        {
            stack.pop_back();
            if(1 < stack.size())
            {
                continue;
            }            
        }
        stack.push_back(*point_it);
        ++point_it;
    }    
    return  stack;  
}
/////////////////////////////////////////////////////////////////////////////////////
T_coord  get_max_triangle_square(const T_points&  points)
{
    T_points  convex_hull_points = get_convex_hull_points(points);
    if(convex_hull_points.size() < 3)
    {
        return 0;
    }
    //Перебираем все тройки выпуклой оболочки, и выбираем с наибольшей площадью треугольника.
    T_coord  S_max = 0;
    for(T_points::const_iterator  it_A = convex_hull_points.begin();
        it_A != convex_hull_points.end(); ++it_A)
    {
        for(T_points::const_iterator  it_B = it_A + 1;
            it_B != convex_hull_points.end(); ++it_B)
        {
            for(T_points::const_iterator  it_C = it_B + 1;
                it_C != convex_hull_points.end(); ++it_C)
            {
                T_point AB = vect(*it_A, *it_B);
                T_point AC = vect(*it_A, *it_C);
                T_coord S = abs(det(AB, AC)) / 2.0;
                if(S_max == 0
                   || S_max < S)
                {
                    S_max = S;
                }
            }        
        }    
    }
    return  S_max;
}
/////////////////////////////////////////////////////////////////////////////////////
int main()
{
    std::locale::global(std::locale(""));
    srand(static_cast<unsigned>(time(0)));
 
    const int N_min = 1;
    const int N_max = 650;
    
    int n = 0;
    do
    {
        std::cout << "n = ";    
        std::cin >> n;
    }while(n < N_min 
           || N_max < n);
    
    T_points  points(n);
    get_points(points);
    std::cout << "Максимальная площадь треугольника из трех точек множества равна "
              << (n < 3 ? 0 : get_max_triangle_square(points))
              << std::endl;
}