Voronoj

#define _CRT_SECURE_NO_WARNINGS
#include <iostream>
#include <assert.h>
#include <math.h>
#include <vector>
#include <algorithm>
#include <queue>

using std::vector;
using std::priority_queue;
using std::cin;
using std::cout;

struct Point {
    double X;
    double Y;
    int Pos;

    Point() : X( 0 ), Y( 0 ), Pos( -1 ) {}
    bool operator==( const Point& other ) const { return X == other.X && Y == other.Y; }
};

struct Vector {
    double X;
    double Y;

    Vector( const Point& p1, const Point& p2 ) : X( p2.X - p1.X ), Y( p2.Y - p1.Y ) {}
    // Знак векторного произведения. Если больше 0, true, то поворот налево.
    bool operator^( const Vector& other ) { return X * other.Y - Y * other.X > 0; }
};

struct Event {
    double Time;
    bool IsNewPoint;
    Point NewPoint;
    int UpperRay;
    int LowerRay;

    Event( const Point& point ) :
        Time( point.X ), IsNewPoint( true ), NewPoint( point ), UpperRay( -1 ), LowerRay( -1 )
    {
    }
    Event( double time, int upperRay, int lowerRay ) :
        Time( time ), IsNewPoint( false ), UpperRay( upperRay ), LowerRay( lowerRay )
    {
    }
    bool operator<( const Event& other ) const { return Time > other.Time; }
};

// Луч, по которому движется пересечение парабол.
struct Ray {
    Point First; // С меньшей x-координатой.
    Point Second; // С большей x-координатой.
    bool IsUp; // true, если над второй точкой.
    bool IsClosed; // Закрыт ли.

    Ray( const Point& first, const Point& second, bool isUp, bool isClosed = false ) :
        First( first ), Second( second ), IsUp( isUp ), IsClosed( isClosed )
    {
    }
    // Y-координата пересечения парабол в момент времени time.
    double GetYAtTime( double time );
};

// Линия - два луча. Если оба луча замкнуты, значит, получился отрезок.
struct Line {
    int Ray1;
    int Ray2;

    Line( int ray1, int ray2 ) : Ray1( ray1 ), Ray2( ray2 ) {}
};

// Дуга параболы.
struct Arc {
    Point Center;
    int UpperRay; // Верхний луч. Индекс в массиве rais. -1, если нет.
    int LowerRay; // Нижний луч. -1, если нет.

    Arc( const Point& point ) : Center( point ), UpperRay( -1 ), LowerRay( -1 ) {}
};

// --------------------------------------------------------------------------------------------------------

double Ray::GetYAtTime( double time )
{
    assert( First.X != Second.X );
    // Решаем систему уравнений:
    // y^2 - 2Fy * y + 2(time - Fx) * x + Fx^2 + Fy^2 - time^2 = 0;
    // x = (Sy - Fy)/(Fx - Sx) * y + (Fx^2 + Fy^2 - Sx^2 - Sy^2)/2(Fx - Sx)
    double xCoef = 2.0 * ( time - First.X );
    double b = -2.0 * First.Y + xCoef * ( Second.Y - First.Y ) / ( First.X - Second.X );
    double c = First.X * First.X + First.Y * First.Y - time * time +
        xCoef * ( First.X * First.X + First.Y * First.Y - Second.X * Second.X - Second.Y * Second.Y ) /
        ( 2.0 * ( First.X - Second.X ) );
    double d = b * b - 4 * c;
    assert( d > 0 );
    d = sqrt( d );
    return IsUp ? ( -b + d ) / 2.0 : ( -b - d ) / 2.0;
}

void ReadPoints( vector<Point>& points )
{
    while( true ) {
        Point p;
        cin >> p.X;
        cin >> p.Y;
        p.Pos = static_cast< int >( points.size() );
        if( cin.eof() ) {
            break;
        }
        points.push_back( p );
    }
}

// -----------------------------------------------------------------------

class CVoronoj {
public:
    CVoronoj( const vector<Point>& _points ) : points( _points ), lines( points.size(), vector<Line>() ) {}

    double Process();

private:
    const vector<Point>& points;
    priority_queue<Event> events;
    vector<Arc> arcs;
    vector<Ray> rais;
    vector<vector<Line>> lines;

    void nextEvent();
    void onNewPoint( double time, const Point& newPoint );
    void onCross( double time, int upperRay, int lowerRay );
    void checkCrossEvent( const Arc& arc );

    double calcAvgPolygonSize();
};

double CVoronoj::Process()
{
    // Инициализируем массив событий.
    for( int i = 0; i < static_cast< int >( points.size() ); ++i ) {
        Event event( points[i] );
        events.push( event );
    }
    // Обрабатываем события.
    while( !events.empty() ) {
        nextEvent();
    }

    return calcAvgPolygonSize();
}

void CVoronoj::nextEvent( )
{
    Event event = events.top( );
    events.pop( );
    if( arcs.empty( ) ) {
        assert( event.IsNewPoint );
        arcs.push_back( Arc( event.NewPoint ) );
        return;
    }
    if( event.IsNewPoint ) {
        onNewPoint( event.Time, event.NewPoint );
    } else {
        // Пересечение.
        onCross( event.Time, event.UpperRay, event.LowerRay );
    }
}

void CVoronoj::onNewPoint( double time, const Point& newPoint )
{
    // Найдем дугу, на которую попадает точка.
    int i = 0;
    while( true ) {
        assert( i < static_cast< int >( arcs.size() ) );
        if( arcs[i].LowerRay == -1 ) {
            break;
        }
        if( rais[arcs[i].LowerRay].GetYAtTime( time ) < newPoint.Y ) {
            // Нашли дугу, содержащую точку.
            break;
        }
        ++i;
    }
    // Добавим два луча одной прямой. И саму прямую.
    Ray upperRay( arcs[i].Center, newPoint, true );
    Ray lowerRay( arcs[i].Center, newPoint, false );
    rais.push_back( upperRay );
    rais.push_back( lowerRay );
    Line line( rais.size() - 2, rais.size() - 1 );
    lines[arcs[i].Center.Pos].push_back( line );
    lines[newPoint.Pos].push_back( line );

    // Разделим дугу arcs[i] на три.
    Arc upper( arcs[i].Center );
    Arc middle( newPoint );
    Arc lower( arcs[i].Center );
    upper.UpperRay = arcs[i].UpperRay;
    upper.LowerRay = rais.size() - 2;
    middle.UpperRay = upper.LowerRay;
    middle.LowerRay = rais.size() - 1;
    lower.UpperRay = middle.LowerRay;
    lower.LowerRay = arcs[i].LowerRay;

    // У верхней и нижней дуги могут появиться пересечения.
    checkCrossEvent( upper );
    checkCrossEvent( lower );

    // Подменим дугу arcs[i] на три в массиве дуг.
    arcs[i] = upper;
    arcs.insert( arcs.begin() + i + 1, 2, middle );
    arcs[i + 2] = lower;
}

void CVoronoj::onCross( double time, int upperRay, int lowerRay )
{
    assert( upperRay != -1 );
    assert( lowerRay != -1 );
    // Пересечение могло запоздать.
    if( rais[upperRay].IsClosed || rais[lowerRay].IsClosed ) {
        return;
    }
    // Найдем дугу, которая схлопнулась.
    int i = 0;
    for( ; i < static_cast< int >( arcs.size() ); ++i ) {
        if( arcs[i].UpperRay == upperRay ) {
            break;
        }
    }
    // Это не может быть первая или последняя дуга.
    assert( i > 0 );
    assert( i < static_cast< int >( arcs.size() - 1 ) );
    assert( arcs[i].LowerRay == lowerRay );

    // Завершим лучи.
    rais[upperRay].IsClosed = true;
    rais[lowerRay].IsClosed = true;

    // Добавим новые лучи.
    Point middleCenter = arcs[i].Center;
    Point upperCenter = rais[upperRay].First == middleCenter ? rais[upperRay].Second : rais[upperRay].First;
    Point lowerCenter = rais[lowerRay].First == middleCenter ? rais[lowerRay].Second : rais[lowerRay].First;
    if( upperCenter.X < lowerCenter.X ) {
        rais.push_back( Ray( upperCenter, lowerCenter, !rais[lowerRay].IsUp, true ) );
        rais.push_back( Ray( upperCenter, lowerCenter, rais[lowerRay].IsUp ) );
    } else {
        rais.push_back( Ray( lowerCenter, upperCenter, !rais[upperRay].IsUp, true ) );
        rais.push_back( Ray( lowerCenter, upperCenter, rais[upperRay].IsUp ) );
    }
    // Добавим прямые.
    Line line( rais.size() - 2, rais.size() - 1 );
    lines[upperCenter.Pos].push_back( line );
    lines[lowerCenter.Pos].push_back( line );

    // Обновим дугу сверху и снизу.
    arcs[i - 1].LowerRay = static_cast< int >( rais.size() ) - 1;
    arcs[i + 1].UpperRay = static_cast< int >( rais.size() ) - 1;
    // Удалим дугу.
    arcs.erase( arcs.begin() + i );
    checkCrossEvent( arcs[i - 1] );
    checkCrossEvent( arcs[i] );
}

void CVoronoj::checkCrossEvent( const Arc& arc )
{
    if( arc.UpperRay == -1 || arc.LowerRay == -1 ) {
        return;
    }
    const Ray& upperRay = rais[arc.UpperRay];
    const Ray& lowerRay = rais[arc.LowerRay];
    Point p2 = arc.Center;
    Point p1 = upperRay.First == p2 ? upperRay.Second : upperRay.First;
    Point p3 = lowerRay.First == p2 ? lowerRay.Second : lowerRay.First;

    // Если точки p1, p2, p3 образуют левый поворот, то пересечение есть.
    if( Vector( p1, p2 ) ^ Vector( p2, p3 ) ) {
        // Найдем время, при котором произойдет пересечение.
        // x = (p2y - p1y)/(p1x - p2x) y + (p1x^2 + p1y^2 - p2x^2 - p2y^2)/2*(p1x - p2x).
        double a1 = ( p2.Y - p1.Y ) / ( p1.X - p2.X );
        double a2 = ( p3.Y - p2.Y ) / ( p2.X - p3.X );
        double b1 = ( p1.X * p1.X + p1.Y * p1.Y - p2.X * p2.X - p2.Y * p2.Y ) / ( 2.0 * ( p1.X - p2.X ) );
        double b2 = ( p2.X * p2.X + p2.Y * p2.Y - p3.X * p3.X - p3.Y * p3.Y ) / ( 2.0 * ( p2.X - p3.X ) );
        double y = ( b2 - b1 ) / ( a1 - a2 );
        double x = a1 * y + b1;
        double dist = sqrt( ( p1.X - x ) * ( p1.X - x ) + ( p1.Y - y ) * ( p1.Y - y ) );
        double time = x + dist;

        Event event( time, arc.UpperRay, arc.LowerRay );
        events.push( event );
    }
}

int main()
{
    vector<Point> points;
    ReadPoints( points );
    CVoronoj voronoj( points );
    cout << voronoj.Process();
    return 0;
}