Разбить вып. многоугольник на 4 равновеликие части (C++), run-time error

[ch-el]

Здравствуйте!
Решаю следующую задачу: на стандартный ввод подаётся число вершин выпуклого многоугольника, а потом пары чисел - координаты вершин в порядке обхода. Необходимо "разрезать" этот многоугольник взаимно перпендикулярными прямыми на 4 равновеликие части - точнее, указать точку, где располагается пересечение этих самых прямых и угол наклона одной из прямой.
Я использую следующий алгоритм, использующий следующие соображения:
1) для любого направления (направлением будем называть пару чисел (A,B), семейство прямых с уравнением Ax+By+C = 0) есть ровно 1 прямая данного направления, делящая многоугольник на 2 равновеликие части (в силу непрерывности);
2) если мы разрежем многоугольник 2 перпендикулярными прямыми, каждая из которых отвечает п.1, то если мы будем перечислять площади кусков по часовой стрелке, то получится последовательность вида {a,b,a,b}, т.е. площади противоположных кусков будут равны, причём такой разрез единственен в силу единственности каждой из прямых.

Произведём разрез 2 прямыми (назовём их первой и второй), каждая из которых делит многоугольник пополам, причём пусть "первая" параллельна оси Ox, и назовём площади получившихся кусков в порядке обхода по часовой стрелке {a, b, a, b}. Не умаляя общности, будем считать a > b. Теперь будем поворачивать прямые по часовой стрелке. Когда угол достигнет pi/2, "первая" прямая будет параллельна уже оси Oy, а последовательность площадей выглядит теперь как {b,a,b,a}, т.е. a и b поменялись ролями. Т.е. первая площадь этой последовательности непрерывно из a перешла в b, т.е. существует такой угол между 0 и pi/2, т.ч. a' = b' => все 4 куска равновелики.

Я написала код на C++, который по идее решает эту задачу таким методом. Если чуть подробнее, то я рассматриваю углы от 0 до pi/2 и методом бисекции ищу нужный - рассматриваю среднее значение угла на текущем отрезке, для этого значения строю прямые по пункту 1, смотрю на разность площадей четвертей, после этого решаю, какую из половин текущего отрезка рассматривать на следующей итерации.

Программа проходит 3 теста с верным ответом, однако на 4ом происходит run-time error, который мне не удаётся, но содержимое теста, к сожалению, неизвестно.
Насколько я понимаю, такая ошибка может возникать, если я повторно очищаю динамически выделенную память (такого вроде не наблюдается), если я выхожу за границы массива (вроде тоже проверила) или делю на 0 (таких возможностей в коде я усмотрела лишь 2 и поставила в тех местах if-ы, отправила код снова и тем не менее run-time error сохранилась). Какие варианты появления этой ошибки ещё возможны?
Кроме того, я нагенерила 1000 выпуклых многоугольников и прогнала код на них - никаких падений программы у меня, по крайней мере, не происходит, хотя я не могу судить, насколько верный ответ она получает.
Ниже привожу код, буду благодарна за помощь в поимке этого рантайма.

[ch-el]

#include <iostream>
#include <cmath>
#include <stdio.h>
#include <stdlib.h>
#include <iomanip>
#include <float.h>


struct Point {
    double xx, yy;
};

int NN;
Point *b_pol;
Point *half_pol;
double b_S;

const double EPS = 1e-7;

double area(Point *pol, int quan);
void find_bisection(int *place, Point *inter, Point line);
double S_difference(double angle, Point& center);

int main()
{
    std::cin >> NN;
    b_pol = new Point[NN + 1];
    half_pol = new Point[NN + 4];
    for (int i = 0; i < NN; ++i) {
        std::cin >> b_pol[i].xx >> b_pol[i].yy;
    }
    b_pol[NN] = b_pol[0];
    b_S = area(b_pol, NN);
    Point center;
    double min_angle = 0, max_angle = M_PI / 2;
    double cur_angle, cur_dif;
    double max_dif = S_difference(max_angle, center);
    double min_dif = S_difference(min_angle, center);
    bool flag = false;
    if (fabs(min_dif) < EPS) {
        flag = true;
        cur_angle = min_angle;
    }
    while (!flag && fabs(min_dif - max_dif) > EPS * b_S) {

        cur_angle = (min_angle + max_angle) / 2.0;
        cur_dif = S_difference(cur_angle, center);
        if (fabs(cur_dif) < EPS) {
            flag = true;
        } else if ((min_dif > DBL_EPSILON && cur_dif < -DBL_EPSILON) ||
                   (min_dif < -DBL_EPSILON && cur_dif > DBL_EPSILON)) {
            max_dif = cur_dif;
            max_angle = cur_angle;
        } else {
            min_dif = cur_dif;
            min_angle = cur_angle;
        }
    }
    printf("%lf %lf\n", center.xx, center.yy);
    printf("%lf\n", (cur_angle * 180) / M_PI);

    delete[] b_pol;
    delete[] half_pol;
    return 0;
}

double area(Point *pol, int quan)
{
    double ans = 0;
    for (int i = 0; i < quan; ++i) {
            ans = ans + (pol[i + 1].xx - pol[i].xx) * (pol[i + 1].yy + pol[i].yy);
    }
    ans = fabs(ans) / 2.0;
    return ans;
}

void find_bisection(int *place, Point *inter, Point line)
{
    int min_ind = 0, max_ind = 0;
    double min_d = line.xx * b_pol[min_ind].xx + line.yy * b_pol[min_ind].yy;
    double max_d = min_d;
    double cur_d;
    for (int i = 0; i < NN; ++i) {
        cur_d = line.xx * b_pol[i].xx + line.yy * b_pol[i].yy;
        if (cur_d > max_d) {
            max_ind = i;
            max_d = cur_d;
        } else if (cur_d < min_d) {
            min_ind = i;
            min_d = cur_d;
        }
    }

    double S_min = 0;
    double S_max = b_S;
    int count = 0;
    double cur_S;
    while (fabs(S_min - S_max) > EPS * b_S) {

        cur_d = (min_d + max_d) / 2.0;
        double test_cur, test_next;
        count = 0;
        for (int side = 0; side < NN; ++side) {
            test_cur = line.xx * b_pol[side].xx + line.yy * b_pol[side].yy - cur_d;
            test_next = line.xx * b_pol[side + 1].xx + line.yy * b_pol[side + 1].yy - cur_d;

            if (test_cur * test_next < -DBL_EPSILON) {
                place[count] = side;
                double A_side = b_pol[side].yy - b_pol[side + 1].yy;
                double B_side = b_pol[side + 1].xx - b_pol[side].xx;
                double C_side = b_pol[side].xx * b_pol[side + 1].yy -
                                b_pol[side + 1].xx * b_pol[side].yy;
                double down = line.xx * B_side - A_side * line.yy;
                if (fabs(down) < EPS) {
                    std::cout << "zero1";
                } else {
                    inter[count].xx = (cur_d * B_side + C_side * line.yy) / down;
                    inter[count].yy = (-line.xx * C_side - A_side * cur_d) / down;
                }
                ++count;
            } else if (fabs(test_cur) < DBL_EPSILON && fabs(test_next) > EPS) {
                place[count] = side;
                inter[count] = b_pol[side];
                ++count;
            }
        }
// A_1 * (B2 - B1) + A2
        half_pol[0] = inter[0];
        count = 1;
        for (int i = place[0] + 1; i <= place[1]; ++i, ++count) {
            half_pol[count] = b_pol[i];
        }
        half_pol[count] = inter[1];
        ++count;
        half_pol[count] = inter[0];
        S_min = area(half_pol, count);
        S_max = b_S - S_min;
        if (min_ind > place[0] && min_ind <= place[1]) {
            if (S_min - S_max > DBL_EPSILON) {
                max_d = cur_d;
            } else {
                min_d = cur_d;
            }
        } else {
            if (S_min - S_max > DBL_EPSILON) {
                min_d = cur_d;
            } else {
                max_d = cur_d;
            }
        }
    }
}

[ch-el]

double S_difference(double angle, Point& center)
{
    int count;
    Point fir, sec;
    fir.xx = sin(angle);
    fir.yy = -cos(angle);
    sec.xx = cos(angle);
    sec.yy = sin(angle);
    Point inter_fir[2], inter_sec[2];
    int place_fir[2], place_sec[2];
    find_bisection(place_fir, inter_fir, fir);
    find_bisection(place_sec, inter_sec, sec);

    double A_fir = inter_fir[0].yy - inter_fir[1].yy;
    double B_fir = inter_fir[1].xx - inter_fir[0].xx;
    double C_fir = inter_fir[0].xx * inter_fir[1].yy - inter_fir[1].xx * inter_fir[0].yy;

    double A_sec = inter_sec[0].yy - inter_sec[1].yy;
    double B_sec = inter_sec[1].xx - inter_sec[0].xx;
    double C_sec = inter_sec[0].xx * inter_sec[1].yy - inter_sec[1].xx * inter_sec[0].yy;

    double down = A_sec * B_fir - A_fir * B_sec;
    if (fabs(down) < EPS) {
        std::cout << "zero2";
    } else {
        center.xx = (C_fir * B_sec - C_sec * B_fir) / down;
        center.yy = (A_fir * C_sec - A_sec * C_fir) / down;
    }

    int stop = 0;
    if (place_fir[0] > place_sec[0]) {
        stop = 1;
    } else if (place_fir[0] == place_sec[0]) {
        double d_fir = (inter_fir[0].xx - b_pol[place_fir[0]].xx) *
                       (inter_fir[0].xx - b_pol[place_fir[0]].xx) +
                       (inter_fir[0].yy - b_pol[place_fir[0]].yy) *
                       (inter_fir[0].yy - b_pol[place_fir[0]].yy);
        double d_sec = (inter_sec[0].xx - b_pol[place_sec[0]].xx) *
                       (inter_sec[0].xx - b_pol[place_sec[0]].xx) +
                       (inter_sec[0].yy - b_pol[place_sec[0]].yy) *
                       (inter_sec[0].yy - b_pol[place_sec[0]].yy);
        if (d_fir - d_sec > DBL_EPSILON) {
            stop = 1;
        }
    }
    half_pol[0] = center;
    half_pol[1] = inter_fir[0];
    count = 2;
    for (int i = place_fir[0] + 1; i <= place_sec[stop]; ++i, ++count) {
        half_pol[count] = b_pol[i];
    }
    half_pol[count] = inter_sec[stop];
    ++count;
    half_pol[count] = center;

    double S_fir = area(half_pol, count);

    half_pol[0] = center;
    half_pol[1] = inter_sec[stop];
    count = 2;
    for (int i = place_sec[stop] + 1; i <= place_fir[1]; ++i, ++count) {
        half_pol[count] = b_pol[i];
    }
    half_pol[count] = inter_fir[1];
    ++count;
    half_pol[count] = center;

    double S_sec = area(half_pol, count);
    return S_fir - S_sec;
}